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The  image  understanding  projects  reported  herein  are  maturing  with 
symbolic  matching,  structure  location,  edge  fitting,  stochastic  texture 
analysis  and  SVD  feature  selection,  all  being  reported  upon  in  some 
detail.  The  image  processing  projects  present  both  new  and  concluding 
projects.  New  projects  include  double  phase  binary  computer  generated 
holograms  and  turntable  radar  imaging  via  coherent  multi-frequency  radar 
return  processing.  Older  projects  resulting  in  successful  theoretics 
and  experimental  work  include  a posteriori  restoration  and  perceptual 
model  color  image  coding.  Our  on-going  smart  sensor  project  is  expandinj 
rapidly  with  old  circuits  being  driven  at  near  real  time  TV  rates  and 
new  circuits  being  designed  for  7x7  area  processing  for  both  enhance- 
ment and  texture  development.  The  Institute  has  installed  a real  time 
TV  solid  state  refresh  monitor  and  display  at  ARPA  headquarters.  This 
allows  recent  pictorial  results  to  be  made  available  over  the  ARPANET. 
Any  and  all  contractors  can  make  use  of  this  device  with  software 
drivers  available  from  the  Institute. 
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Abstract 


This  technical  report  summarizes  the  image 
understanding,  smart  sensor,  and  image  processing  research 
activities  performed  by  the  Image  Processing  Institute  at 
the  University  of  Southern  California  during  the  period  of  1 
October  1977  through  31  March  1978  under  Contract  Number 
F-33615-76-C-1 203  with  the  Advanced  Research  Projects  Agency 
Information  Processing  Techniques  Office. 

The  image  understanding  projects  reported  herein  are 
maturing  with  symbolic  matching,  structure  location,  edge 
fitting,  stochastic  texture  analysis  and  SVD  feature 
selection,  all  being  reported  upon  in  some  detail.  The 
image  processing  projects  present  both  new  and  concluding 
projects.  New  projects  include  double  phase  binary  computer 
generated  holograms  and  turntable  radar  imaging  via  coherent 
multi-frequency  radar  return  processing.  Older  projects 
resulting  in  successful  theoretic  and  experimental  work 
include  a posteriori  restoration  and  perceptual  model  color 
image  coding.  Our  on-going  smart  sensor  project  is 
expanding  rapidly  with  old  circuits  being  driven  at  near 
real  time  TV  rates  and  new  circuits  being  designed  for  7x7 
area  processing  for  both  enhancement  and  texture 
development.  The  Institute  has  installed  a real  time  TV 
solid  state  refresh  monitor  and  display  at  ARPA 
headquarters.  This  allows  recent  pictorial  results  to  be 
made  available  over  the  ARPANET.  Any  end  all  contractors 
can  make  use  of  this  device  with  software  drivers  available 
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Research  Overview 


The  past  six  months  have  been  quite  productive  on  a 
variety  of  research  and  development  fronts.  The  image 
understanding  projects  are  maturing  with  symbolic  matching, 
structure  location,  edge  fitting,  stochastic  texture 
analysis  and  SVD  feature  selection,  all  being  reported  upon 
in  some  detail.  The  image  processing  projects  present  both 
new  and  concluding  projects.  New  projects  include  double 
phase  binary  computer  generated  holograms  and  turntable 
radar  imaging  via  coherent  multi-frequency  radar  return 
processing.  Older  projects  resulting  in  successful 
theoretic  and  experimental  work  include  a posteriori 
restoration  and  perceptual  model  color  image  coding.  Our 
on-going  smart  sensor  project  is  expanding  rapidly  with  old 
circuits  being  driven  at  near  real  time  TV  rates  and  new 
circuits  being  designed  for  7x7  area  processing  for  both 
enhancement  and  texture  development.  The  Institute  has 
recently  acquired  a high  precision  hardcopy  color  device  for 
improved  output  capability  and  has  installed  a real  time  TV 
solid  state  refresh  monitor  and  display  at  ARPA 
headquarters.  This  allows  recent  pictorial  results  to  be 
made  available  over  the  ARPANET.  Any  and  all  contractors 
can  make  use  of  this  device  with  software  drivers  available 
from  the  Institute.  Finally  this  past  six  months  has 
witnessed  the  graduation  of  one  Ph.D.  student  and  numerous 
Institute  personnel  publications  listed  at  the  end  of  this 
semi-annual  report. 


2 ■ Image  Understanding  Projects 


The  projects  reported  herein  describe  progress  on  a 
variety  of  fronts  in  our  image  understanding  efforts. 
Segment  matching  is  discussed  in  some  detail  followed  by 
symbolic  matching  of  images  with  substantial  changes  in 
orientation.  Results  in  both  these  areas  are  encouraging  as 
the  reader  will  see  in  reading  these  contributions.  A 
technique  is  next  presented  for  locating  structures  in 
aerial  images  followed  by  the  presentation  of  a new  edge 
fitting  algorithm.  A quantitative  analysis  tool  is 
developed  to  more  accurately  define  what  and  where  an  edge 
is.  These  contributions  are  followed  by  a discussion  on 
stochastic  texture  analysis  in  which  some  previously 
accepted  perceptual  moment  assumptions  are  challenged. 
Finally  the  section  closes  with  a discussion  on  singular 
value  decomposition  image  feature  extraction  as  a useful 
tool  for  obtaining  features  in  an  image  understanding 
scenario. 


2.1  Matching  Segments  of  Images 
Keith  E.  Price 


Introduction 

Image  registration  has  few  direct  applications  by 
itself,  but  it  is  a necessary  step  for  computer  change 
analysis,  tracking  with  multiple  images,  or  motion  analysis. 
We  define  registration  as:  finding  corresponding  points  or 
regions  in  two  or  more  images  of  a single  scene,  where  the 
images  may  be  taken  from  a different  position  or  at  a 
different  time.  For  some  registration  and  change  analysis 
systems,  one  image  may  be  transformed  so  that  each  point  in 
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one  image  precisely  corresponds  with  the  same  point  in  the 
other  image,  but  this  is  not  always  necessary. 

There  are  two  basic  types  of  registration  methods: 
signal  based  and  symbolic.  Signal  based  registration 
systems  find  pairs  of  corresponding  points  by  matching  the 
image  values  in  a small  neighborhood  around  the  two  points 
using,  for  example,  a correlation  function  [1-6].  A 
symbolic  registration  system  locates  pairs  of  corresponding 
segments  by  comparing  values  of  features  over  entire 
segments  [7-10] . The  symbolic  representation  of  an  image 
contains  segments  representing  objects  or  part  of  objects 
and  a feature  based  description  of  the  segments. 

Signal  based  registration  systems  must  make  a major 
assumption  about  the  allowable  pairs  of  images:  the  two 
images  are  already  in  close  alignment  and  major  portions  of 
the  scene  remain  unchanged.  The  first  is  necessary  to  limit 
the  area  of  the  image  where  corresponding  points  may  be 
located  to  reduce  the  search  time,  and  the  second  assures 
that  enough  corresponding  points  for  a good  match  will  look 
the  same  in  both  images.  The  first  restriction  can  be 
relaxed  if  the  approximate  camera  positions  are  known  so 
that  one  image  can  be  transformed  before  any  other 
processing  is  attempted.  These  assumptions  are  valid  for 
controlled  stereo  pairs  and  thus  most  stereo  analysis  can  be 
performed  very  well  at  the  signal  level.  Image  based 
matching  is  not  always  possible  for  arbitrarily  selected 
points  if  the  point  falls  in  the  middle  of  a homogeneous 
region,  thus  some  processing  to  locate  edges  or  other 
features  of  regions  would  be  necessary  in  a general  system 
[6]  . 


Analysis  after  registration  and  some  change  analysis  is 
usually  best  performed  using  symbolic  techniques,  e.g. 
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shadow  detection  cannot  be  performed  on  a single  pixel,  but 
requires  the  analysis  of  nearby  pixels  to  determine  if  there 
might  be  a shadow  and  the  location  of  an  object  that  could 
cast  a shadow  in  that  position.  Thus  it  would  seem 
practical  to  attempt  to  perform  the  initial  registration 
symbolically  so  that  the  change  analysis  and  later 
processing  can  more  easily  be  performed  symbolically. 

Related  work  has  been  reported  in  [8]  and  [9]  where  a 
limited  feature  based  description  is  used  to  compare  images 
of  a moving  scene  and  maintain  a description  of  the  scene. 
The  description  is  limited  to  features  invariant  to 
rotations  and  translations,  and  is  also  limited  by  the  fact 
that  the  scenes  are  composed  of  white  objects  on  a black 
background.  This  system  uses  a model  of  the  objects  in 
motion  in  the  scene  to  predict  the  current  view  so  that  it 
can  match  scenes  with  occlusions.  The  assumptions  made  and 
problems  studied  in  [9]  differ  from  those  in  the  work 
reported  here  so  the  results  are  not  directly  comparable, 
but  both  represent  similar  techniques  being  applied  to 
different  aspects  of  the  image  matching  problem. 

We  will  first  describe  the  problem  which  we  are 
attempting  to  solve  and  present  the  scenes  which  will  be 
used  to  illustrate  important  aspects  of  the  solution.  The 
third  section  will  describe  the  procedures  which  we  have 
implemented  to  solve  these  problems.  Section  IV  will 
present  some  of  the  results  from  these  programs. 

Problem  Description 

We  are  exploring  the  general  problem  of  finding  changes 
between  different  views  of  the  same  scene.  The  views  may  be 
taken  at  different  times  or  from  different  positions,  and 
the  changes  can  occur  in  any  of  features  used  to  describe 
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objects  in  the  scene  - size,  shape,  color,  location,  etc. 

As  a necessary  step  in  locating  these  differences  we  have 
worked  on  a system  to  locate  corresponding  segments  in  two 
views  of  one  scene.  The  differences,  in  many  cases,  are 
then  simply  the  differences  in  the  feature  values  of  the  two 
corresponding  segments.  Several  assumptions  have  been  made 
about  the  type  of  input  and  the  extent  of  changes  which 
occur.  First,  the  images  are  arbitrary  views  of  natural 
scenes,  including  ground  level  views  and  aerial  photographs, 
and  have  already  been  segmented  by  a computer  program  such 

as  in  (7,  11-13].  This  implies  that  the  segmentations  are 

not  perfect  and  may  vary  from  one  view  to  the  next.  Second, 
objects  may  appear  in  only  one  view,  and  many  objects  may  be 
very  similar  to  each  other.  Also,  there  may  be  both  global 
and  local  changes  in  the  scene  and  changes  can  occur  in  any 
feature,  but  most  corresponding  region  pairs  will  be  more 
similar  than  noncorresponding  region  pairs.  Finally, 
included  in  the  descriptions  of  the  tasks  given  to  this 

program  is  an  indication  of  which  features  will  usually 

remain  constant  and  which  features  will  usually  change.  The 
uses  of  these  assumptions  should  become  apparent  in  the  next 
section . 

Scene  Descriptions 

Figures  1 through  4 show  the  initial  images  and  partial 
segmentations  for  an  urban  scene.  There  are  size  and 
position  changes  for  the  objects  between  the  two  images  but 
no  rotation  differences.  There  are  many  local  differences 
between  the  images  and  these  are  discussed  further  in  [7], 

In  this  study  we  are  interested  only  in  locating 
corresponding  positions  in  the  two  images  and  determining 
the  scale  difference.  This  task  can  be  accomplished  best  by 
using  several  well  defined,  easily  segmented  regions, 
therefore  only  the  bright  regions  have  been  segmented  in 
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Figure  3.  Urban  scene,  segmentation  of  Image  1 

The  segmented  regions  are  the  gray  areas 
in  the  figure. 


Figure  4.  Urban  scene,  segmentation  of  Image  2 
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this  pair  of  images. 

Figures  5 through  8 give  the  input  images  and 
segmentations  for  a rural  scene.  These  images  are  slightly 
rotated  with  respect  to  each  other  but  were  taken  at  about 
the  same  time.  The  two  segmentations  include  only  the  large 
untextured  (homogeneous)  regions  and  the  much  smaller  bright 
regions.  This  scene  illustrates  the  operation  of  the 
matching  procedure  when  there  is  a small  orientation 
difference  and  when  there  are  several  similar  regions 
extracted  in  the  segmentation  operation.  These  two  are  also 
discussed  further  in  [7]  along  with  the  description  of  the 
operation  of  the  segmentation  procedure  as  applied  to  both 
of  these  scenes. 

Matching 

In  a system  which  is  designed  to  match  a pair  of  real 
images,  the  problem  is  usually  not  whether  the  segments 
match  exactly,  using  a given  feature,  but  how  close  the 
segments  match.  When  applied  to  the  problem  of  finding  the 
corresponding  region  in  a set  of  regions,  the  question 
becomes:  how  well  do  these  two  regions  match  compared  to  how 
well  the  region  in  the  first  image  matches  other  regions  in 
the  second  image?  Thus  we  want  a match  function  which  will 
indicate  how  well  two  segments  match  by  combining  how  well 
segments  match  for  each  feature. 

A first  approximation  of  how  well  two  regions  match  for 
one  feature  is  given  by  the  difference  between  the  feature 
values.  This  might  be  reasonable  if  all  features  had 
similar  values  (e.g.  only  spectral  features),  but  when 
combining  arbitrary  features  the  absolute  difference  is  not 
meaningful.  Therefore,  the  feature  value  differences  must 
be  weighted  so  that  the  contribution  to  the  matching 
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Figure  7 


1 


Rural  scene,  segmentation  of  Image  1 
Each  of  the  large  regions  is  displayed  with  a 
different  gray  level  and  the  small  bright 
regions  are  given  one  gray  value. 


Figure  8.  Rural  scene,  segmentation  of  Image  2 


function  of  each  feature  is  approximately  the  same  for 
equally  poor  matches.  This  feature  normalization  value  must 
account  for  several  variations:  the  sizes  of  the  feature 
values  are  very  different  - some  range  from  0 to  1,  some 
from  10  to  one  million  (size) ; and  some  features  vary  more 
with  small  differences  in  the  image  than  others  or  are  more 
prone  to  errors  in  computation.  Table  1 gives  a list  of  the 
features  and  the  weighting  for  the  actual  matching 
procedure.  We  chose  a simple  method  to  combine  the  weighted 
differences,  just  the  sum  of  the  weighted  absolute 
differences.  Thus  the  region  to  region  match  rating  is: 


n 


Rij  - - E ivik-V\ 

k=l 


(1) 


where  R. . is  the  value  of  the  match  between  regions  i and  j, 
ij 

n is  the  number  of  valid  features,  V^^V.^)  is  the  value  of 
the  kth  feature  of  regions  i(j),  and  W^  is  the  normalization 
weight  for  the  kth  feature.  Region  i and  j are  from 
different  images.  The  result  is  negative  so  that  better 
matches  have  a greater  value,  i.e.  closer  to  zero. 


The  neighbor  and  relative  position  features  have  no 
affect  until  at  least  one  pair  of  corresponding  regions  is 
located.  The  neighbor  feature  value  for  region  i and  j are 
computed  as  follows:  Let  be  the  set  of  all  regions  which 
are  neighbors  of  region  i and  already  have  a corresponding 
region  in  the  second  image.  Let  Nj  be  the  set  of  the 
corresponding  regions  for  those  regions  in  N^ . Then  let  Nj 
be  the  members  of  Nj  which  are  neighbors  of  region  j.  The 
neighbor  feature  value  for  region  i is  the  length  of  and 

for  region  j is  the  length  of  N, . The  relative  position 

/ t J 

feature  (above,  below,  to  right,  and  to  left)  value  is 
defined  in  an  equivalent  manner. 
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Fea  tur e 


Feature  Weighting 


Size 

Location,  each  coordinate 
Per imeter 2/Area 
Or ienta  tion 
Length  to  Width  Ratio 
Area/Area  of  Minimum 
Bounding  Rectangle 
Area/Area  of  Bounding 
Ellipse 

Color  or  Intensity 

Neighbors 

Relative  Positions 


Minimum  of:  5/Area  of  region, 

100/Area  of  image 
8/Image  dimension 
2A/P2 
2 
2 

1.6667 

2.5 

l/2o  Standard  deviation  of  the  color 
parameter  over  the  giver,  region. 
2/Number  of  neighbors  with  a 
corresponding  region 
2/Number  of  relative  regions  with  a 
corresponding  region 


Table  1.  FEATURE  VALUE  NORMALIZATION  WEIGHTS.  These  are  the  weights 
used  in  Eqs.  (1)  and  (2)  for  all  the  available  features.  These 
values  are  used  to  make  the  contribution  of  each  feature  to  the  match 
rating  approximately  equal.  When  feature  values  are  indicated  in  the 
definition  of  the  weighting,  the  values  are  from  image  1 or  the  region 
in  image  1.  The  features  with  constant  weight  have  values  in  a 
limited  range,  such  as  0 to  1. 
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We  stated  in  the  problem  description  section  that  the 
task  description  may  include  the  fact  that  some  features  are 
changing  and  are  thus  unreliable  for  use  in  the  matching 
process.  Alternatively,  we  may  be  given  the  fact  that  some 
feature  must  not  change  between  the  two  views  and  these 
features  are  very  reliable.  Therefore,  for  a given  task, 
some  features  should  contribute  more  to  the  region  to  region 
match  rating  than  other  features  and  some  features  may  not 
even  be  considered.  We  added  a strength  parameter  to  the 
rating  function  to  account  for  the  different  types  of 
features.  This  gives  a region  to  region  match  rating  of: 

n 

- - Z |Vlk-Vjk|MkSk 
k=l 

where  is  the  strength  of  the  kth  feature.  There  are  only 
three  possible  strength  values  whicn  are  used  by  the 
program.  These  values  were  selected  so  that  it  requires 
several  poor  matches  in  lower  strength  features  to  have  the 
same  impact  as  an  equally  poor  match  in  one  higher  strength 
feature.  The  features  are  assigned  to  strengths  based  on 
information  given  in  the  task  description,  with  features 
assigned  the  medium  strength  when  the  outside  information 
does  not  specify  otherwise.  If  a good  match  in  a feature  is 
indicated  as  absolutely  necessary  then  the  feature  is  given 
the  higher  strength  and  if  it  is  given  that  the  feature  may 
change,  but  not  too  dramatically,  then  the  lower  strength  is 
used.  Features  which  are  not  avail  abl  e and  those  which  are 
known  to  be  totally  unreliable  for  this  task  are  not 
considered  in  the  rating  - effectively  a fourth  strength  of 
zero.  Generally,  the  assignment  of  strengths  is  the  same 
for  all  matches  between  a pair  of  images. 

Certain  global  changes  may  mean  that  some  features 
cannot  be  used  in  the  matching  operation,  for  example  scale 
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changes  mean  that  size,  and  absolute  position  will  not 
match.  But,  adjustments  in  the  feature  values  can  be  made 
if  these  global  changes  are  known.  We  indicated  earlier 
that  knowledge  of  which  properties  are  changing  may  be 
given,  but  not  the  magnitude  of  these  changes,  which  would 
be  necessary  to  make  any  adjustments.  But  detailed 
information  about  the  global  changes  is  available,  if 
several  corresponding  region  pairs  can  be  located  without 
using  these  changing  features.  After  the  adjustments  are 
computed  the  feature  can  be  used  like  any  other  feature, 
except  that  the  feature  values  will  be  corrected  before  they 
are  used  in  the  matching  function.  In  this  case,  the 
strength  assignments  will  change  after  the  first  few  matches 
because  the  adjusted  feature  values  will  match  better  than 
the  original  feature  values. 

This  region  matching  process  is  the  basic  operation  for 
the  procedure  which  locates  corresponding  regions  in  two 
images.  The  symbolic  registration  procedure  compares  the 
given  region  with  each  potential  matching  region  using  given 
sets  of  features  with  different  strengths  (e.g.  those  which 
must  remain  constant  and  those  which  may  change) . The 
region  which  best  matches  the  given  region  is  returned  as 
the  corresponding  region.  It  should  be  noted  that  this  is 
not  an  attempt  to  find  a complete  one  to  one  mapping  from 
regions  in  one  image  to  regions  in  another  image,  since  this 
is  not  possible  when  objects  may  appear  or  disappear  and 
different  regions  may  be  segmented  in  the  two  images.  The 
regions  to  be  registered  and  the  strengths  assigned  to  each 
feature  are  not  selected  by  the  symbolic  registration 
procedure,  but  are  given  by  other  elements  of  the  total 
system  (7).  The  order  in  which  the  corresponding  regions 
are  located  may  be  important.  It  is  especially  important 
that  initial  matches  be  accurate  if  the  neighbor  and 
relative  position  features  are  to  be  of  any  value.  The 
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usual  order  is  to  find  matches  for  the  regions  in  order  of 
size  so  that  the  neighbor  and  relative  position  features  can 
be  used.  But  when  a few  initial  matches  are  needed  to  make 
adjustments  for  global  changes,  the  large  regions  may 
provide  unreliable  information  expecially  if  they  are 
adjacent  to  the  edge  of  the  image.  Because  of  this,  easily 
identified  regions,  other  than  the  larger  regions,  may  also 
be  matched  first,  but  these  decisions  are  made  a controlling 
system  or  the  user,  and  are  based  on  the  task  description. 

Results 

Figures  9 and  10  show  the  registration  results  for  the 
two  scenes,  i.e.  the  corresponding  regions  in  the  two 
images  - the  corresponding  regions  have  the  same  intensity 
values  in  the  two  pictures  which  make  up  each  figure.  Since 
only  the  bright  regions  were  segmented  in  the  urban  scene, 
these  are  the  only  ones  which  are  used  in  the  registration 
process.  This  pair  of  images  has  a global  scale  and 
position  change  which  is  computed  after  the  first  few 
corresponding  regions  are  located  (the  large  horizontal 
object  in  the  lower  right  and  the  left  most  round  object  in 
the  upper  center).  These  adjustments  are  used  for  all 
subsequent  matching  operations  so  that  location  and  size  can 
also  be  treated  as  constant  features.  This  adjustment  is 
necessary  to  obtain  correct  matches  for  many  of  the  regions, 
especially  the  group  of  round  regions  near  the  top  of  the 
scene  where  the  shape  of  the  extracted  segments  differ 
because  of  shadows  which  make  position  an  important  features 
in  the  matching  operation.  The  computed  global  size  change 
is  about  1.5. 

The  large  regions  in  the  rural  scene  images  are  large 
homogeneous  (nontextur ed)  regions  and  there  are  many 
differences  in  the  segmentations  of  these  regions  because 
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Figure  9.  Symbolic  Registration  for  urban  scene 

The  corresponding  regions  have  the  same  gray 
value  in  the  two  pictures.  Image  1 is  on  the 
left  and  image  2 on  the  right.  Some  of  the 
small  round  regions  may  not  be  visible  since 
their  values  are  almost  white. 


Figure  10.  Symbolic  Registration  for  rural  scene 

Corresponding  regions  have  the  same  gray  level 
in  the  two  pictures.  Image  1 is  on  the  left 
and  image  2 on  the  right.  The  bright  areas 
near  the  center  are  almost  white  and  some  may 
not  be  apparent  in  the  pictures. 
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this  feature  is  very  sensitive  to  noise.  Because  of  this, 
there  are  some  regions  which  have  no  correct  corresponding 
region  among  the  segmented  regions  in  the  other  image.  But 
generally,  when  they  exist,  the  correct  corresponding 
regions  are  located.  The  results  for  the  smaller  regions 
might  be  improved  if  an  adjustment  for  the  rotation 
difference  is  calculated  and  used,  but  many  of  the  houses 
are  correctly  matched  even  though  they  are  all  similar  to 
each  other. 

One  immediate  question  that  may  be  raised  is  how  close 
the  best  match  rating  is  to  ratings  for  other  potential 
matches.  The  mean  rating  for  the  best  match  for  the  rural 
scene  shown  here  is  -287  (a=116)  and  the  mean  for  the  second 
best  matching  region  (i.e.  the  region  which  would  have  been 
selected  if  the  actual  matching  region  had  not  been 
segmented)  is  -530  (a=301) . Values  similar  to  these  are 
common  for  other  scenes  which  we  have  analyzed.  In  most 
cases  the  mean  for  the  second  best  matching  pairs  is  about 
one  standard  deviation  from  the  mean  for  the  best  matching 
pairs  (in  some  cases  it  is  much  further,  but  rarely  much 
closer).  Tables  2 and  3 show  more  detailed  region  to  region 
matching  results  for  the  urban  and  rural  scene  respectively. 
Regions  in  the  first  image  correspond  to  rows  in  the  table, 
and  those  in  image  2 to  columns.  The  values  are  the  region 
to  region  match  values  for  generated  by  comparing  regions  in 
image  1 to  regions  in  image  2.  The  values  along  the 
diagonal  are  the  match  values  for  the  computed  corresponding 
regions.  The  regions  in  each  table  are  ordered  by  the  size 
of  the  region  in  the  first  image.  There  is  a wide  disparity 
between  the  values  of  accepted  matches,  from  -77  to  -590. 
The  differences  between  the  rating  for  the  accepted  pair  and 
other  pairs  can  be  very  small  when  the  regions  are  almost 
identical  such  as  the  round  regions  in  the  urban  scene, 
noted  by  *,  and  the  houses  in  the  rural  scene,  marked  by  +. 
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-77  -1369  -925  -1251  -791  -1285  -1166  -1156  -1200  -1201  -1107  -1281  -1182  -1125 

-1735  -410  -1436  -741  -1359  -1061  -569  -925  -1244  -1139  -1252  -1258  -1156  -1249 

-1002  -1150  -300  -1075  -366  -1416  -994  -1182  -1712  -1517  -1670  -1768  -1606  -1574 

-1500  -561  -1126  -277  -1189  -984  -467  -808  -1235  -1060  -1203  -1211  -1108  -1184 

-874  -1120  -430  -1101  -119  -1251  -867  -989  -1549  -1350  -1505  -1590  -1445  -1414 
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To  reduce  the  size  of  this  table  all  numbers  have  been  divided  by  10  and  the  minus  signs  have  been 
removed,  otherwise  the  format  is  the  same  as  Table  2.  The  " + "s  indicate  the  regions  which  correspond 
to  houses  in  the  scene. 


In  some  cases  two  regions  in  the  first  image  match  with  the 
same  region  in  the  second  image  so  some  columns  are  equal. 

Conclusions 

This  work  is  intended  to  be  part  of  a general  purpose 
analysis  system  and  not  simply  the  solution  to  a particular 
set  of  problems.  The  symbolic  registration  procedure  worked 
well  with  the  machine  segmentations  and  performance  should 
improve  if  segmentations  were  improved.  This  research  shows 
that  symbolic  techniques  can  be  used  for  analysis  of  a 
larger  class  of  scenes  than  is  possible  with  signal  based 
system,  but  may  be  used  in  conjunction  with  signal  based 
systems  where  detailed  analysis  of  changes  in  individual 
pixels  is  required.  Since  many  change  analysis  tasks  are 
best  performed  at  a symbolic  level,  it  seems  reasonable  to 
perform  the  matching  and  registration  symbolically  rather 
than  generate  symbolic  descriptions  of  unreliable  difference 
images,  thus  symbolic  registration  will  be  an  important  step 
of  general  change  analysis  systems. 
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2.2  Symbolic  Matching  and  Analysis  with  Substantial  Changes 
in  Orientation 

Keith  E.  Price 


Introduction 


Several  methods  have  been  developed  which  can  be  used 
to  find  correspondences  in  pairs  of  images  of  a changing 
scene  [1-61 . But,  for  various  reasons,  these  systems  did 
not  operate  on  images  which  had  major  changes  in 
orientation  - unless  an  approximate  value  for  the  difference 
was  known  a priori.  Image  based  methods  [1-3]  were  not  used 
to  attempt  a solution  to  this  problem  because  of  the 
complexity  of  searching  for  corresponding  points.  Other 
work  [5,6]  assumes  that  the  views  are  close  - in  time  and 
position  - so  that  major  changes  in  the  point  of  view  are 
not  possible.  But  in  a more  general  image  matching  and 
analysis  system  [4J,  this  type  of  problem  must  be 
considered . 

We  have  undertaken  a series  of  experiments  to  examine 
whether  the  matching  techniques  described  in  [4]  can  be 
easily  applied  to  pairs  of  images  which  have  substantial 
global  differences.  In  this  paper,  we  present  the  results 
of  applying  this  basic  technique  on  pairs  of  images  with 
orientation  differences  of  45°,  90°  and  180°.  These 
orientation  changes  are  in  addition  to  other,  less  drastic. 
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changes  which  may  occur  between  the  two  views.  The  point  is 
to  show  that  symbolic  techniques  can  be  applied  when  there 
are  substantial  global  changes  in  addition  to  the  normal 
local  changes  which  occur  between  two  views  of  one  scene. 

We  will  first  discuss  the  images  which  will  be  used  for 
this  experiment  and  describe  the  results  which  are  desired. 
Then  we  will  present  an  outline  of  the  symbolic  matching 
procedure  [4]  and  some  initial  results  using  this  method. 
The  results  suggest  some  modifications  which  are  then 
described,  followed  by  more  extensive  results  using  the 
modified  procedure. 

Tasks  for  Matching 

The  pairs  of  images  which  we  will  use  here  have  been 
used  earlier  [7],  but  in  the  previous  analyses  there  were  no 
substantial  global  changes.  We  start  with  a pair  of  images 
of  a scene  taken  from  slightly  different  positions,  and 
generate  the  orientation  changes  by  rotating  the  digital 
representations  of  the  second  image  of  the  pair.  The 
original  first  image  and  the  rotated  second  images  are  then 
processed  to  generate  symbolic  descriptions  (a  segmentation 
into  distinct  parts  plus  a description  of  the  segments) 
which  are  used  by  the  matching  procedures.  The  details  of 
the  segmentation  and  description  are  given  elsewhere  and  are 
not  important  for  this  paper  [4,7,8]. 

The  basic  task  to  be  executed  for  the  images  presented 
here  is  to  find  the  corresponding  regions  in  the  two  images 
when  there  is  a large,  but  unknown,  orientation  difference. 
A by-product  of  this  matching  should  be  some  indication  of 
the  actual  orientation  difference.  The  images  in  this 
experiment  were  rotated  45°,  90°,  and  180°,  but  the  features 
which  are  orientation  dependent  are  not  used  even  if  they 
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are  independent  to  rotations  of  90°  (e.g.  ratio  of  area  and 
area  of  minimum  bounding  rectangle) . 

The  first  pair  of  images,  a house,  is  snown  in 
Figure  1.  These  are  color  images,  shown  here  in  black  and 
white,  so  there  are  several  spectral  features  available  for 
use  in  the  segmentation,  description,  and  matching 
operations.  Figure  2 gives  the  segmentations  of  the 
original  first  view  and  the  three  rotated  second  views. 
There  are  some  differences  between  the  regions  segmented  in 
the  first  view  and  each  of  the  second  views,  and  a few 
differences  among  the  three  second  views.  The  images  are 
not  square  so  the  display  program  puts  a white  band  on  the 
right  or  bottom , depending  on  which  dimension  is  smaller. 
Additionally  unsegmented  areas  are  displayed  as  white  areas. 

The  second  pair  of  images  is  shown  in  Figure  3.  These 
are  side  looking  radar  views  and  there  is  only  one  spectral 
input,  so  that  a good  segmentation  is  more  difficult  and  the 
description  is  less  detailed  than  the  preceding  house  scene. 
The  important  changes  in  this  scene  are  in  objects  which  are 
too  small  to  be  segmented  by  our  general  segmentation 
techniques  but  might  be  easily  located  by  special  methods. 
Figure  4 gives  the  four  segmentations  which  will  be  used. 
There  are  a few  large  regions  segmented  in  the  first  view 
which  are  not  segmented  in  the  second  views.  These 
differently  textured  regions  vary  in  size  and  appearance 
between  the  two  views  and  would  not  be  used  in  the  matching, 
therefore  they  were  not  segmented  in  the  second  views.  The 
regions  which  are  segmented  are  the  untextured  areas,  all 
uniformly  dark,  some  of  which  appear  the  same  in  both  views. 
The  locations  of  these  few  regions  ca/i  be  used  to  determine 
the  global  changes  and  could  aid  another  system  in  locating 
the  detailed  changes. 
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Matching  Procedure  Outline 


A more  detailed  description  of  the  matching  procedure 
has  been  presented  elsewhere  [4]  and  will  only  be  outlined 
here.  This  matching  procedure  uses  a feature  based, 
symbolic,  description  of  the  images.  The  basic  unit  in  this 
description  is  an  individual  segment  generated  by  an 
automatic  segmentation  procedure  [8].  These  segments  are 
usually  regions  in  the  images,  but  linear  features  can  be 
described,  too.  Features  which  characterize  properties  such 
as  color,  texture,  size,  shape,  position,  and  adjacencies  are 
used . 

The  matching  procedure  is  also  given  an  indication  of 
which  features  are  available  for  matching  the  current  pair 
of  images,  and  what  strength  to  give  to  the  mismatch  using 
each  of  these  feature.  For  example,  some  features  are  not 
always  computed,  red  and  green  in  a black  and  white  picture, 
and  some  are  given  as  more  likely  to  change  than  others  and 
are  thus  given  less  weight  in  the  matching  operation.  The 
match  procedure  computes  a rating  for  the  match  between  two 
regions  which  is  essentially  a sum  of  the  weighted  absolute 
differences  of  feature  values.  The  weights  used  in  the  sum 
are  composed  of  a normalization  factor  to  make  the 
contribution  from  each  feature  approximately  equal  and  a 
strength  factor  to  account  for  the  different  strengths 
assigned  each  feature.  There  is  only  a small  set  of 
possible  strength  values,  currently  3 different  values. 

Known  global  changes  between  the  two  images  can  be  used 
to  adjust  some  feature  values,  such  as  size,  position  and 
orientation.  But,  these  changes  are  not  given  a priori  and 
must  be  computed  from  a few  initial  pairs  of  matching 
regions.  Thus,  the  very  clearly  defined  regions  should  be 
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Fig.  5.  Initial  matching  regions  for  house  first  view  (left) 
and  house  second  view  rotated  90°  (right) . 


matched  first  so  that  they  may  be  used  for  calculating  some 
global  changes.  In  this  context,  clearly  defined  regions 
are  regions  with  extreme  values  for  some  feature,  e.g. 
largest,  brightest,  longest,  etc. 

We  will  now  present  results  of  applying  this  procedure 
to  the  rotated  images  and  discuss  what  changes  are  necessary 
to  achieve  accurate  results. 

Initial  Results 


Figures  5 and  6 show  the  results  obtained  with  ’the 
above  matching  procedure  for  the  house  scene  (90°  and  180°). 
In  these,  and  all  other  figures,  the  corresponding  regions 
are  displayed  at  the  same  intensity  in  the  pair  of  output 
pictures.  Similar  results  are  obtained  for  45°  and  for 
matching  the  second  image  to  the  first,  but  the  point  here 
is  to  point  out  some  of  the  problems.  The  orientation 
feature  adjustment  was  computed,  using  the  sky  or  roof,  and 
was  also  used  to  get  these  results,  but  there  are  still  many 
errors.  Most  of  the  unmatched  regions  would  match  to  an 
incorrect  corresponding  region  if  a match  is  attempted. 

The  major  problem  is  that  the  location  of  the  region  is 
needed  to  correctly  locate  matches  for  many  of  the  smaller 
regions.  This  is  especially  the  case  when  there  are  size 
and  shape  changes  due  to  segmentation  differences,  such  as 
in  the  bush,  window  and  door  regions.  If  exact  camera 
transformations  are  known  then  the  locations  in  one  image 
can  be  mapped  exactly  into  locations  in  the  other,  but  this 
transformation  is  not  known.  We  are  using  many  features 
other  than  position  so  an  appropriate  mapping  is  sufficient 
to  allow  the  use  of  the  absolute  position  features.  Given  3 
pairs  of  corresponding  regions  we  can  compute  a 
transformation  which  will  map  coordinates  in  one  image  to 


Fig.  5.  Initial  matching  regions  for  house  first  view  (left) 
and  house  second  view  rotated  90°  (right) . 
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Fig.  6.  Initial  matching  regions  for  house  first  view  (left) 
and  house  second  view  rotated  180°  (right) . 
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coordinates  in  the  other  by  solving  2 sets  of  3 equations 
and  3 unknowns.  This  transformation  is  not  optimal  for  all 
regions  in  the  image  and  only  accounts  for  rigid,  global 
changes  - e.g.  rotations  and  translations.  But  this 
transformation  does  make  the  position  features  useable  when 
there  are  large  global  orientation  changes. 

Final  Results 

Figures  7-9  show  the  results  for  the  house  matching 

using  the  computed  location  transformations  - a different 

transformation  is  computed  for  each  image  pair.  Results  are 

given  for  matching  regions  in  image  1 with  image  2 and 

regions  in  image  2 with  image  1.  In  this  set  of  images  the 

sky,  lawn,  and  one  of  the  wall  sections  are  the  initial  3 

regions  used  for  the  transformation  computation.  The 

transformations  do  not  rotate  the  coordinates  precisely  45°, 

90°,  or  180°  because  of  differences  in  segmentations  (the 

sky  and  lawn  are  adjacent  to  the  edge  and  this  cause  some 

changes  in  size  and  shape)  and  a small  orientation 

difference  which  existed  before  the  large  rotations  were 

added.  But,  the  adjustment  is  accurate  enough  to  use  the 

o 

location  feature  in  the  match  operation.  The  results  at  45 
are  less  accurate  than  for  the  other  two,  but  most  of  the 
extra  mistakes  are  accounted  for  by  the  greater  differences 
in  segmentations  (see  Fig.  2).  When  two  regions  in  one 
image  correspond  to  one  region  in  the  second  image,  such  as 
the  2 "bushes"  on  the  right  side  of  the  house  in  the  second 
image  appearing  as  one  region  in  the  first  image,  only  one 
correspondence  appears  in  the  output. 


Figures  10  and  11  present  the  results  for  the  radar 
images.  There  are  very  few  (i.e.  6)  corresponding  regions 
in  these  two  images  and  two  of  the  pairs  have  very  large 
size  changes.  The  results  for  1 8 cP  are  identical  to  those 
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(b)  House  2 rotated  45°  with  house  1. 


Fig.  7.  Matching. 
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(b)  House  2 rotated  90°  with  house  1. 


Fig.  8.  Hatching 
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(b)  House  2 rotated  180°  with  house  1. 


Fig.  9.  Matching 
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(b)  Radar  2 rotated  45°  with  radar  1. 


Fig.  10.  Matching 
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for  90°  and  are  not  presented.  Two  of  the  corresponding 
pairs  for  the  image  1 to  image  2 matching  may  be  difficult 
to  see  since  they  are  nearly  white  in  the  results  picture 

one  is  a correct  match  (the  reversed  "C"  shape  in  the  lower 

left) , and  the  other  is  incorrect  (a  blob  near  the  top 
right).  The  reverse  "C"  region,  the  river  (lower  right)  and 
the  blob  above  the  river  are  the  three  regions  used  to 
compute  the  transformation  in  this  set  of  images.  This 
scene  shows  that  when  these  symbolic  techniques  are  applied 
to  scenes  with  a reduced  feature  set  - no  colors  and  no 

neighboring  regions,  accurate  results  are  possible. 

Summary 

The  complete  symbolic  registration  system  presented 
here  has  the  following  basic  steps: 

1.  Segment  both  images  of  the  scene. 

2.  Generate  a feature  based  description  of  the 
segmented  images. 

3.  Find  corresponding  regions  for  the  most  obvious 
regions. 

4.  Set  orientation  and  size  correction  factors,  if 
necessary . 

5.  Find  several  corresponding  region  pairs. 

6.  Compute  an  approximate  coordinate  transformation, 
if  necessary. 

7.  Using  transformed  positions,  find  all  corresponding 
region  pairs. 

The  matching  results  depend  somewhat  on  the  quality  of 
the  segmentation,  but  the  results  of  these  experiments  show 
that  this  symbolic  technique  can  be  used  to  find 
corresponding  areas  in  pairs  of  images  even  when  there  are 
major  global  changes.  We  would  expect  similar,  or  better, 
results  for  pairs  of  images  with  global  scale,  position,  and 
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scenes 


color  changes.  We  expect  less  reliable  results  for 
with  major  global  changes  in  all  four  (orientation,  scale, 
color,  and  position)  because  so  few  features  are  invariant 
to  all  these  changes  (e.g.  relative  size,  shape  measures, 
and  neighbors).  But  if  a controlling  system  could  provide 
proper  guidance,  corresponding  regions  might  be  located 
which  would  account  for  each  of  the  global  changes, 
separately.  For  example,  scale  changes  could  be  based  on 
matching  the  largest  regions,  orientation  changes  might  be 
based  on  regions  with  distinctive  shape,  and  so  on.  But, 
primary  regions  with  unusual,  or  extreme,  feature  values 
could  be  used  when  there  are  many  global  changes.  In 
conclusion,  symbolic  matching  methods  can  work  with  major 
global  differences,  these  differences  can  be  detected,  and 
they  can  be  used  to  great  advantage  in  later  analysis. 
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2.3  Locating  Structures  in  Aerial  Images 
Ramakant  Nevatia  and  Keith  E.  Price 


(Some  of  this  work  was  reported  previously  in  the  last 
progress  report  and  presented  at  the  October,  1977  ARPA 
workshop.  This  material  is  repeated  here  for  completeness 
and  new  results  are  reported.) 

This  section  describes  our  initial  attempts  to  utilize 
our  segmentation  techniques  to  perform  higher  level  tasks  of 
locating  user  specified  structures  in  aerial  images.  We  use 
edge  based  [1]  as  well  as  region  based  segmentation 
techniques  [2-3],  as  experience  indicate  that  the  two 
methods  have  complementary  strengths  [4].  The  edge 
techniques  are  suitable  for  extracting  elongated  features 
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1 


such  as  roads  and  the  region  methods  for  homogeneous  areas 
such  as  lakes. 

Analysis  of  aerial  images  is,  in  general,  a complex 
task.  The  reasons  for  such  complexities  are  many  and 
varied.  A prime  cause  is  the  presence  of  texture  which 
causes  difficulties  for  the  low  level  processes  such  as  edge 
detection  and  segmentation.  Another  source  of  difficulty  is 
that  the  desired  objects  and  structures  may  be  small 
compared  to  the  size  of  a complete  image.  A detailed 
analysis  of  a complete  high  resolution  aerial  image  is 
generally  prohibitive  because  of  the  computational  costs. 

For  many  applications,  however,  a complete  and  general 
analysis  is  unnecessary.  Specific  structures  of  interest 
may  have  special  properties,  known  a priori,  that  allow  for 
their  easy  extraction.  The  problem  of  searching  for  small 
structures  is  helped  by  locating  them  by  their  spatial 
relationships  to  larger,  more  easily  located  structures.  We 
assume  that  the  user  specifies  the  properties  useful  for 
locating  the  desired  structures  and  also  the  properties  of 
and  relations  with  other  objects.  An  interactive  system  to 
facilitate  the  entry  of  such  data  is  also  being  developed. 
Such  a priori  information  is  likely  to  be  available  in  many 
applications  where  different  images  of  the  same  area  need  to 
be  analyzed. 

Model  and  Descr iption 

The  a priori  information  is  represented  as  a graph 
structure  with  the  objects  represented  as  nodes  and 
relationships  between  objects  as  arcs.  Generally,  the 
properties  of  objects  and  relationships  between  objects 
should  be  unrestricted.  Currently  an  object  is  described  as 
a collection  of  related  line  segments  and  regions.  Both 
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lines  and  regions  are  described  by  position  (both  absolute 

and  relative),  size,  and  orientation  features,  and  regions 

2 

also  have  shape  (parameter  /area,  length/width,  area/area  of 
minimum  bounding  rectangle),  color  (intensity,  red,  green, 
etc.),  and  texture  (e.g.  variance)  properties.  The 
description  and  matching  operations  use  both  regions  and 
lines  with  no  direct  indication  of  the  type  of  segment 
needed,  because  the  lines  are  essentially  described  as 
narrow  elongated  regions. 

Several  relationships  between  regions  are  currently 
used  such  as  neighbors,  nearby  regions,  and  the  symbolic 
relationships  of  above,  below,  left  and  right.  Other 
relationships  such  as  symmetry  and  similarity  are  useful  but 
have  not  been  implemented.  Simpler  relations  such  as 
larger,  brighter,  etc.  are  not  directly  implemented,  but 
could  be  easily  computed  from  given  feature  values. 

The  image  is  described  in  the  same  manner  as  the  user 
specified  objects,  except  that  the  image  description  is 
performed  automatically,  without  direct  user  input.  We  will 
use  the  term  model  to  mean  the  descriptions  of  the  user 
defined  object  or  grqup  of  objects  which  we  are  attempting 
to  locate  in  an  image. 

For  example,  in  Figure  1 we  have  an  aerial  view  of  a 
mostly  rural  scene  and  we  wish  to  locate  the  docks  along  the 
lower  side  of  the  river.  It  might  be  possible  to  search  the 
entire  image  looking  for  the  desired  objects,  but  it  seems 
to  be  easier  to  locate  the  docks  by  using  larger  regions  and 
linear  features  such  as  lakes  and  roads.  Figure  2 shows  a 
hand  drawn  pictorial  representation  of  this  portion  of  the 
image  with  only  the  important  objects  represented.  Figure  3 
gives  the  internal  graph  structure  representation  of  objects 
which  can  be  used  to  locate  the  docks.  The  a priori 


-43- 


knowledge  must  also  provide  the  system  with  a way  of 
locating  the  area  containing  the  docks  given  the  positions 
of  the  other  related  regions  - the  river,  lakes,  and  road. 


Our  representation  and  use  of  knowledge  is  similar  to 
that  described  by  Tenenbaum  [5].  The  main  difference  is  in 
Tenenbaum's  use  of  single  pixel  attributes  to  uniquely 
distinguish  objects  in  a given  context.  We  use  the  object 
attributes  to  aid  in  the  segmentation  of  the  image  and  then 
use  the  attributes  of  large  segmented  parts  for  recognition. 
Bajcsy  and  Lieberman  [6]  attempted  to  model  general  outdoor 
scenes,  but  the  detailed  models  which  we  use  are  specific  to 
the  particular  task  being  performed  and  scene  being 
analyzed.  Bajcsy  and  Tavakoli  [7]  used  special  purpose 
programs  to  locate  specific  objects  (i.e.  roads  and 

bridges)  rather  than  a general  matching  technique. 

Matching  Strategy 

For  simplicity,  and  because  it  proved  to  be  accurate, 
we  use  a general  purpose  symbolic  matching  system  to  find 
the  elements  in  the  image  which  corresDond  to  the  individual 
elements  in  the  user  specified  model  [3].  The  general 
matching  procedure  locates  corresponding  regions  by 

comparing  the  symbolic  description  of  a region  in  one  image 
with  the  symbolic  descriptions  of  regions  in  another  image, 
and  selecting  the  region  which  is  most  similar.  The 
similarity  function  for  two  regions  is  essentially  the  sum 
of  normalized  absolute  differences  of  the  feature  values  of 
the  two  regions.  This  gives  a region  to  region  match 
function  of: 

n 

Ru  - 'vik  - Hi 


n is  the  number  of  valid  features,  V^(V^ 
the  kth  feature  of  regions  i(j),  and  is 
weight  for  the  kth  feature.  Region  i 
different  images.  The  result  is  negative 
matches  are  larger,  i.e.  closer  to  zero. 


) is  the  value  of 
the  normalization 
and  j are  from 
so  that  better 


In  the  general  image  matching  case,  the  user  can 
specify  that  some  features  will  change,  or  are  less  reliable 
for  the  given  task.  Those  features  are  then  given  less 
weight  in  the  computation  of  the  similarity  function. 
Conversely  other  features  can  be  given  as  required  matches 
and  these  are  given  more  weight  in  the  rating.  These 
task-dependent  weights  for  the  feature  value  differences  are 
in  addition  to  the  normalization  factors  which  are  constant 
for  all  tasks.  The  weights  were  selected  so  that  a poor 
match  in  a highly  weighted  feature,  i.e.  a constant 
feature,  will  have  the  same  impact  as  several  poor  matches 
in  lower  weighted  features.  The  different  weights  should  be 
less  important  when  matching  with  the  user  specified  model, 
since  the  given  features  are  only  the  ones  important  for  the 
description  and  thus  are  features  which  should  not  change. 

The  system  locates  the  best  corresponding  segment  in 
the  image  for  each  segment  in  the  model,  but  no  attempt  is 
made  to  find  the  globally  best  one-to-one  mapping.  Because 
of  this,  the  order  in  which  corresponding  segments  are 
located  may  be  important.  Since  we  are  using  features  such 
as  neighboring  regions,  and  regions  above,  below,  left  and 
right,  an  incorrect  match  early  in  the  process  will  lead  to 
inaccurate  similarity  ratings  later.  The  default  order  for 
matching  is  by  size,  but  other  orders  can  be  specified  by 
the  user  or  made  available  automatically  if  desired.  We 
assume  that  the  ordering  information  is  supplied  by  the 
user,  and  have  made  no  attempt  to  automate  the  strategy 
generation  process  as  in  [8]. 
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We  now  present  results  of  applying  these  procedures  on 
two  scenes. 

Results 


The  first  scene  is  the  one  presented  in  Figures  1 to  3 
in  the  previous  section.  This  is  an  aerial  view  of  a mostly 
rural  area.  The  goal,  for  this  example,  is  to  locate  the 
docks  along  the  lower  edge  of  the  large  river.  Figure  2 
gives  a pictorial  description  of  the  important  elements  in 
the  scene,  such  as  might  appear  on  a map  or  a sketched 
description  of  the  desired  objects.  The  portion  of  the 
image  where  the  docks  can  be  located  is  specified  in 
relation  to  the  objects  in  the  model  description.  Thus  the 
bounds  of  this  area  can  be  computed  when  the  position  in  the 
image  of  the  objects  in  the  model  are  known.  This  limited 
area  can  be  further  analyzed  to  complete  the  task.  Figure  4 
shows  the  lines  and  regions  which  are  segmented  in  this 
image,  because  the  lakes  and  river  are  the  only  regions  used 
in  the  description.  These  are  the  only  ones  extracted  by 
the  segmentation  operation. 

When  the  user  specified  model  is  matched  with  the 
machine  segmented  image,  the  regions  and  line  shown  in 
Figure  5 are  selected  as  the  corresponding  segments  in  the 
scene.  The  location  and  dimensions  of  the  upper  lake  and 
the  location  of  the  lower  edge  of  the  river  are  used  to 
specify  the  area  where  the  docks  are  located.  This  area  is 
extracted  from  the  complex  image  for  detailed  analysis  with 
which  the  individual  docks  may  be  located.  The  two 
subimages  corresponding  to  the  dock  area  are  given  in 
Figure  6. 

The  results  for  this  scene  require  a few  comments.  The 
segmentation  procedures  failed  to  separate  the  lower  lake 
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Figure  4.  Segmentation  of  Image  1. 


Figure  5 


Regions  Corresponding  to  Model  Elements . 


from  the  river  region  since  the  two  regions  are  connected  by 
a narrow  neck,  but  a separation  should  be  possible  with  more 
processing  of  the  extracted  regions  to  eliminate  narrow 
necks.  In  this  case  the  neck  was  just  wide  enough  to  escape 
the  processing  which  is  already  applied.  Because  of  this, 
there  is  no  actual  region  in  the  image  which  corresponds  to 
the  lower  lake  in  the  model.  The  region  in  the  image  which 
is  selected  as  the  corresponding  region  is  the  upper  lake 
region.  In  the  procedure  which  computes  the  position  of  the 
dock  areas,  the  position  of  the  best  matching  lake  is  used. 
Because  of  the  errors  in  the  relative  position  of  the 
selected  lower  lake  region,  i.e.  it  is  on  the  wrong  side  of 
the  road,  the  upper  lake  was  matched  best  and  thus  this  is 
the  region  which  is  used.  This  illustrates  the 
insensitivity  of  the  higher  level  processing  to  certain 
errors  in  the  lower  level  processing. 

We  also  applied  the  same  system  to  another  scene,  shown 
in  Figure  7.  This  is  an  aerial  photograph  dominated  by  a 
river,  a central  power  plant,  a large  group  of  houses  and 
several  prominent  roads.  The  task  for  this  scene  is  to 
locate  the  power  plant  area.  This  image  has  about 
2000  x 2000  points  so  that  a detailed  analysis  of  the  entire 
image  would  be  impractical.  The  position  of  the  power  plant 
is  specified  by  the  location  of  certain  roads  and  river 
segments.  Figure  8 illustrates  the  important  objects  in  the 
scene  which  we  wish  to  locate.  The  location  of  the  river 
regions  and  some  of  the  roads  is  important  mainly  because  of 
the  relation  of  these  segments  to  other,  necessary,  road 
segments.  Figure  9 shows  the  segments  extracted  from  this 
image.  The  segments  include  the  large  untextured  regions, 
small  bright  regions,  and  long  line  features.  Figure  10 
shows  the  regions  selected  by  the  system  as  those  which 
correspond  to  the  elements  in  the  user  defined  model.  The 
position  of  the  centrally  located  road  segment  is  used  to 
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specify  the  subwindow  which  includes  the  power  plant  and 
very  little  else.  Figure  11  shows  the  subwindow  which  was 
extracted  from  this  image.  This  subwindow  is  512  x 512 
pixels  and  could  be  further  analyzed  by  this  system  to  find 
specific  objects  if  a description  of  this  portion  of  the 
scene  is  furnished  by  the  user. 

Other  experiments  were  tried  with  this  image  and  model 
to  determine  which  features  are  necessary  for  locating  the 
regions  in  the  model.  Correct  results  were  obtained  using 
only  the  size  and  relative  positions  (above,  below,  etc.)  of 
the  segments.  This  is  accounted  for  mostly  by  the  fact  that 
the  sizes  of  the  regions  are  very  large  compared  to  other 
regions  (the  lines  are  also  the  longest  lines  which  were 
extracted) , and  the  sizes  specified  in  the  model  were  close 
to  the  sizes  in  the  image.  If  the  sizes  in  the  model  were 
less  precise  then  other  features  would  be  necessary. 

Conclusions 


Some  results  of  processing  a complex  aerial  image  using 
both  the  line  and  the  region  based  techniques  have  been 
shown.  It  appears  that  the  use  of  simple  techniques, 
specifically  suited  to  particular  objects  in  an  image,  may 
allow  useful  processing  of  rather  complex  images.  At  this 
time  the  array  of  segmentation  attributes  is  limited.  While 
it  is  hoped  that  the  described  techniques  have  general 
applicability,  our  experience  with  real  images  is,  as  yet, 
limited. 
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2.4  A Mew  Edge  Fitting  Algorithm 
Ikram  E.  Abdou 


Introduction 

In  a previous  paper  [1]  it  was  mentioned  that  the 
Hueckel  operator  is  not  an  optimum  edge  fitting  method 
because  the  operator  is  based  on  some  approximations  which 
affect  the  final  performance.  This  paper  introduces  a new 
edge  fitting  algorithm  which  avoids  some  of  the 
disadvantages  of  the  Hueckel  operator.  A discussion  of  the 
basic  ideas  of  the  new  algorithm  and  a comparison  between 
its  performance  and  the  Hueckel  performance  are  given  in  the 
following  sections. 

Basic  Concepts 

In  edge  fitting  the  image  function  f(i,j)  defined  over 
a subregion  J)  is  compared  with  an  ideal  edge  model  S (i,j), 

tO 

where  to  is  a set  of  edge  parameters.  The  difference  between 

the  actual  and  ideal  models  is  function  of  to  , and  by 

changing  these  parameters  the  difference  can  be  minimized. 

The  edge  acceptance  is  based  on  the  value  of  the  minimum 

difference.  If  it  is  less  than  a given  threshold  (t)  , the 

image  subregion  is  classified  as  an  edge  with  the 

corresponding  parameters  to  . . Usually  the  mean  square 

min 

error  is  used  to  measure  the  difference  between  the  ideal 
and  actual  edge.  This  error  is  given  in  the  form 

E„  - D E m 

j 1 

- 

The  ideal  edge  model  used  in  the  new  algorithm  is  a 
central  edge  with  slope  along  one  of  the  four  basic 


directions  (shown  in  Figure  lb) . For  this  model  w is  the 
set  {0^,n,a,A}  , where 

0^  is  the  edge  orientation 

n is  the  edge  width 

a is  the  edge  average  height 

A is  the  slope  of  the  central  ramp, 

is  minimized  by  changing  two  discrete  variables,  0^  and 
n,  and  two  continuous  variables,  a and  A.  Since  0^  and  n 
assume  finite  number  of  values,  the  minimization  problem  can 
be  solved  by  repeating  the  computation  for  each  discrete 
variable  and  choosing  the  minimum  E, , for  all  and  n.  As 
an  example,  assume  that  f(i,j)  is  fitted  to  the  vertical 
edge  shown  in  Figure  la,  Sw(i,j)  is  now  given  by 


Su)(i’j)  = a_An 

= a+Ai 
= a+An 

Substituting  in  Eg.  1. 


-N  £ i < -n,  V j 
-n  £ i £ n,  V j 
n < i < N,  V j 


- (n+1) 

n 

E = V 

0)  / ^ 

yt  (a-An-f (i,j))2+  XI  (a+A 

i=-N 

j i=-n  j 

Expanding  and  arranging  terms.  E 

-(n+1) 

- Z E 

(a-f (i, j))2+A  { 2n  Sf(i* 

* j 

i=-N  T 

N 

, -(n+1) 

-2n 

Z «i.J>}+*2{  Z En2+ 

i=n+l 

j -N  j 

= CQ  + CXA 

+ C242 
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i=-n  3 


i=n+l  j (3. a) 


(3.b) 


For  given  edge  orientation  and  width,  E is  minimum 

0) 

when 

8E 

air  = 0 (4*a) 

3E 

y _ n 

3A  U ( 4 . b) 

Substituting  in  Eq.  (3) 


(2N+1) ' 


(5. a) 


i J 


A " 


(5 . b) 


where  C]  and  C2  are  as  given  in  Eq.  3.  In  addition 


E . 


min 


c 


(6) 


Similar  expressions  can  be  obtained  for  different  edge 
orientation  and  widths.  The  actual  image  is  fitted  to  the 
edge  model  which  results  in  minimum  error.  In  the 
evaluation  of  Eq.  6 it  is  not  necessary  to  calculate  Cg 
because  it  is  constant  for  the  same  image  subregion.  The 
acceptance  of  the  edge  can  be  based  on  the  values  of  Amin 
and  E^in • The  computation  time  required  by  the  previous 
procedure  is  much  less  than  the  Hueckel  time.  While  the 
results  obtained  are  more  accurate  as  will  be  shown  in  the 
following  section. 
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Exper lmental  Results 

The  performance  of  the  edge  fitting  algorithm  with 
different  mask  sizes  is  evaluated  using  the  figure  of  merit 
defined  in  (2).  The  results  obtained  are  shown  in  Figure  2. 
Comparing  these  results  with  the  results  obtained  for  the 
Hueckel  operator  (1),  it  is  clear  that  the  new  algorithm  has 
better  performance. 

Another  question  of  interest  is  to  what  extent  is  the 
new  edge  fitting  operator  better  than  a 3-level  simple 
operator  (3)  of  the  same  mask  size?  The  performance  curves 
of  3-level  simple  operator  with  mask  sizes  5,  and  9 are 
shown  in  Figure  2.  It  is  clear  that  for  small  mask  size  and 
very  low  SNR  the  edge  fitting  algorithm  is  not  so  good  as 
simple  masks  operator.  This  observation  can  be  explained  by 
the  fact  that  the  edge  fitting  algorithm  bases  its  decision 
on  an  estimation  of  the  edge  parameters.  This  estimation  is 
sensitive  to  noise  specially  when  the  number  of  pixels  used 
is  small.  However  the  edge  fitting  operator  has  better 
performance  for  high  SNR  and  for  large  mask  sizes. 
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Figure  2.  Edge  location  figure  of  Merit  as  a function 
of  signal-to-noise  ratio. 
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2.5  Stochastic  Texture  Analysis 
William  K.  Pratt* 


In  the  early  Sixties,  Julesz  (1)  attempted  to  determine 
the  parameters  of  stochastic  texture  fields  of  perceptual 
importance.  These  experiments  were  extended  by  Julesz  et 
al.  (2-4)  about  ten  years  later.  There  has  also  been 
recent  work  consisting  of  further  extensions  and  variations 
of  Julesz's  original  experiments  by  Pollack  [5]  and  by  Purks 
and  Richards  [6],  One  difficulty  in  application  of  these 
experiments  to  image  classification  and  analysis  is  that  the 
stochastic  fields  investigated  did  not  adequately  simulate 
natural  texture.  The  fields  were  only  quantized  to  2 to  4 
gray  levels,  and  they  only  possessed  one-dimensional 
correlation.  The  following  summarizes  the  results  of  a 
recent  study  by  Pratt,  Faugeras,  and  Gagalowicz  [7]  on  the 
analysis  of  stochastic  texture  fields  containing  multiple 
gray  levels  and  possessing  two-dimensional  correlation. 

Figure  1 contains  a block  diagram  for  a general  model 
of  stochastic  texture  generation.  An  array  of  independent, 
identically  distributed  pixels  W(j,k)  with  probability 
density  p(W)  passes  through  a linear  or  nonlinear  system 
described  by  the  operation  (>{•}  to  produce  the  stochastic 
texture  field  F(j,k).  Reference  [7]  provides  design  methods 
for  selecting  p(W)  andO{»}  in  order  to  achieve  textural 
fields  with  desired  statistical  properties. 


*This  work  was  performed  in  collaboration  with  O.D. 
Faugeras  and  A.  Gagalowicz  of  Institut  de  Recherche 
d ' Informatique  et  d ' Automatique , Le  Chesnay,  France. 
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Figure  2 contains  computer-generated  texture  fields  for 
an  experiment  to  determine  perceptual  sensitivity  to 
differences  in  probability  density  between  texture  field 
pairs.  Moments  listed  on  the  figure  are  defined  elsewhere. 
Figure  2a  offers  a comparison  of  256  x 256  pixel  fields 
of  256  gray  levels  with  different  independent  first  order 
densities  but  identical  means.  The  A field  possesses  a 
uniform  density,  while  the  B field  has  a density  shaped  like 
an  isosceles  triangle.  In  figure  2b,  there  is  a split  field 
comparison  between  correlated  fields  with  common  uniform 
first  densities,  but  different  second  order  densities. 
Visual  discrimination  is  readily  apparent  in  both  examples. 
Figure  2c  contains  photographs  of  texture  pairs  in  which 
pixels  are  pairwise  independent.  But,  their  third  order 
probability  densities  differ.  Here,  discrimination  is  not 
possible.  Experiments  of  this  nature  have  led  Julesz  to 
conjecture  that  texture  fields  differing  only  in  third  and 
higher  order  densities  cannot  be  discriminated  by  a human 
observer.  One  possible  explanation  for  the  lack  of 
discrimination  in  figure  2c  is  the  absence  of  spatial 
correlation.  However,  this  is  not  the  proper  explanation. 
In  figure  2d  the  random  numbers  of  figure  2c  have  been 
spatially  correlated  by  separable  linear  processing  along 
rows  and  columns  to  produce  fields  with  first  order 
Markovian  correlation,  but  differing  third  order  densities. 
Since  discrimination  is  not  possible,  Julesz's  conjecture 
has  been  confirmed  for  spatially  correlated  as  well  as 
uncorrelated  fields. 

The  experiments  of  figure  3 illustrate  visual 
sensitivity  to  changes  in  the  autocorrelation  of  Gaussian 
density  texture  fields  with  identical  first  order  densities. 
Texture  differences  corresponding  to  changes  of  about  10%  in 
the  correlation  lag  parameters  are  visible. 


-67- 


(c)  Uncorrelated  third 
order  density  fields 

oA  = 0.289  aB  = 0.289 

aA  = 0.000  oig  = 0.000 

6.  = 0.058  eD  = -0.058 


(d)  Correlated  third  order 
density  fields 

o.  = 0.167  = 0.167 


Field  comparison  of  texture  fields  with  differing 
probability  densities.  n*  = nR  = 0.500. 


Figure  4 proves  the  insufficiency  of  the 
autocorrelation  function  as  a description  of  texture.  In 
figures  4a  and  4b  the  fields  are  generated  with  identical 
means,  variances,  and  autocorrelation  functions,  but 
different  third  order  densities.  The  texture  differences 
are  easily  discernible. 

The  experiments  of  figures  2 to  4 serve  to  establish 
useful-  bounds  for  developing  stochastic-based  visual  texture 
features.  Second  order  statistical  measures  should  be 
sufficient,  but  the  mean,  variance,  and  autocorrelation 
function  measures,  by  themselves,  although  directly  or 
indirectly  necessary,  are  not  sufficient.  The  task  then  is 
to  determine  those  stochastic  features  within  the 
established  bounds  that  are  perceptually  sufficient. 

The  stochastic  texture  examples  presented  here  are 
important  not  only  for  the  inherent  information  they  convey 
about  visual  texture  perception;  they  also  provide  valuable 
test  fields  for  evaluation  of  candidate  texture  features.  A 
universally  sufficient  texture  feature  vector  must  provide 
sufficient  information  to  quantitatively  discriminate  the 
Gaussian  texture  fields  of  figure  3 possessing  differing 
autocorrelation  functions,  and  also  discriminate  the 
non-Gaussian  fields  of  figure  4 that  have  identical 
autocorrelation.  If  such  feature  vectors  pass  this 
stringent  test,  then  they  are  worth  evaluating  for  natural 
texture. 
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2.6  Singular  Value  Decomposition  Image  Feature  Extraction 
Behnam  Ashjari  and  William  K.  Pratt 


The  singular  value  decomposition  (SVD)  is  a numerical 
technique  of  matrix  transformation  by  which  an  arbitrary 
matrix  can  be  expressed  in  outer  product  form.  SVD 
<•«:  ir'ii'mi  have  been  applied  in  the  solutions  of 
. « * ; >ne>J  seta  of  linear  equations  (1).  There  has 


also  been  recent  application  of  the  SVD  concept  to  image 
restoration  and  coding  [2,3] . Another  application,  explored 
here,  is  the  use  of  the  SVD  to  extract  descriptive  features 
from  an  image  region  for  purposes  of  classification  or 
analysis . 


Singular  Value  Matrix  Decomposition 


Consider 
non-negative , 
defined  to  be 


an  N x N matrix  F of  rank  R containing 
real  elements.  The  SVD  transform  of  F is 


where  U and  V are  unitary  matrices  and  S is 
matrix  whose  diagonal 

(1)  > A^(2)  > ...  > A^(R)  ■ ...\(N) 


a diagonal 
elements 
are  called  the 

singular  values  of  F.  Since  U and  V are  unitary  matrices, 
the  inverse  transform  yields 

T (2) 

F = US V1 


The  unitary  matrix  U 

contains  columns  u 

-n 

composed  of 

T 

eigenvectors  of  the 

symmetric  matrix 

product 

m 

FFA  . 

Similarly,  the  columns  v 

of  V are  eigenvectors 
n — 

of  r f. 

The 

defining  relations  are 

T 

SS 

= UTtFFT]U  = A 

(3a) 

T 

SS 

= VTfFTF]V  = A 

(3b) 

Thus,  the  squares  of  the  singular  values 
diagonal  matrix  i\  are  the  eigenvalues  of 
Consequently,  the  singular  values  can 
generated.  From  solution  of  the 

characteristic  equation 


arranged  in  the 
both  FTT  and  FTF. 
be  algebraically 
N roots  of  the 
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N 

C ( A ) = E aixj  = detTFF1  - XI] 

j =0  J V ' 

It  should  be  observed  that  the  eigenvalues  _A  are  positive, 

T T ~ 

real  var’ables  since  FF  and  F F are  real,  symmetric 
matrices.  For  consistency,  the  singular  values  S are  chosen 
to  be  the  positive  square  roots  of  the  corresponding 
eigenvalues . 


It  is  often  useful  to  express  the  matrix  decomposition 
of  eq. (2)  in  the  outer  product  expansion 

R 

F = E S(J)  u,vT  (5) 

j = l J J 

where  the  sum  is  over  the  rank  R of  F.  In  a similar  manner, 
the  symmetric  matrix  products  can  be  expressed  in  the  outer 
product  forms 


N 


FFT  = E 

T 

f f1 

(6a) 

m=l 

— m-fn 

m N 

pTp  . ^ 

h hT 
— m— m 

(6b) 

m=l 


where  f denotes  the  m-th  column  of  F and  h is  the  m-th  row 


SVD  Image  Feature  Extraction  Concept 

The  singular  values  of  a matrix  can  be  considered  as 
descriptors  or  features  of  the  matrix  elements  and  their 
inter-relationships.  If  the  matrix  is  composed  of  randomly 
chosen  real  numbers,  the  singular  values  will  tend  toward 
equality.  On  the  other  hand,  a highly  structured  matrix 
will  exhibit  a few  dominating  singular  values.  Figure  1 
sketches  the  qualitative  structural  relationship  of  the 
singular  values.  This  observation  forms  the  basis  for 
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utilization  of  the  SVD  as  a means 
An  SVD  expansion  is  formed  over  a 
some  measure  of  the  skewness 
distribution  is  formed  to 
"coherence"  of  the  pixel  values, 
solidified  shortly. 


of  forming  image  features, 
block  of  N x N pixels,  and 
of  the  singular  value 
characterize  the  spatial 
This  vague  concept  will  be 


Deterministic  Properties 


Certain  deterministic  image  patterns  possess  SVD 


expansions  that  can  be  easily  computed.  For  example,  an 
array  of  unit  value  pixels  can  be  generated  by  the  outer 
product  expansion 


[1  1 • 


• • 1] 


P 1 * 
11. 


11 


11.. 


(7) 


Such  a matrix  is  of  unit  rank,  and  therefore,  possesses  only 
one  non-zero  singular  value.  A matrix  containing 
alternating  vertical  stripes  of  zeroes  can  be  formed  by 


*1* 

[1010.  . .10] 

'10  10.  . .10' 

1 

10  10.  . .10 

....  • • 

= 

. ...  • • 

....  . • 

.1. 

_1  0 1 o . . .10. 

(8) 


Horizontally  striped  arrays  can  be  generated  by  reversing 
the  positions  of  the  vectors  in  eq.(8).  The  reason  that  the 
arrays  of  eqs. (7)  and  (8)  have  only  one  non-zero  singular 
value  is  that  their  rows  are  linearly  dependent,  in  fact 
identical.  Hence,  any  array  with  repeated  rows  (or  columns) 
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will  possess  only  one  non-zero  singular  value.  This  class 
includes  striped  arrays  of  various  periods  and  arrays  with 
single  horizontal  or  vertical  line  strokes. 

A checkerboard  array  of  ones  and  zeroes  can  be 
generated  by  the  sum  of  two  outer  products 

[10  10.  . .10] 

+ 


and  therefore,  possesses  two  non-zero  singular  values. 
Finally,  an  N x N identity  matrix  with  ones  down  its  main 
diagonal  and  zeroes  elsewhere  is  a rank  N matrix  and  has  N 
equal  singular  values. 

These  simple  deterministic  examples  illustrate  some 
interesting  points.  First,  the  SVD  expansion  is  invariant 
to  the  period  of  striped  horizontal  and  vertical  patterns. 
In  fact,  a constant  amplitude  array,  which  can  be  considered 
an  array  with  a 100%  duty  cycle  period,  has  the  same  number 
of  non-zero  singular  values  (one)  as  a striped  matrix  with 
single  pixel  stripes.  As  a consequence,  the  SVD  is  clearly 
not  useful  as  a measure  of  periodic  structure  of 
deterministic  arrays.  The  second  major  point  of  the 
examples  is  that  rotation  of  a pattern  affects  the  number  of 
singular  values.  A checkerboard  pattern  consisting  of 
diagonal  stripes  possesses  two  non-zero  singular  values, 
while  horizontally  or  vertically  striped  arrays  only  have 
one  non-zero  singular  value.  As  an  even  more  extreme 
example  consider  the  N x N diagonal  matrix  with  a single 
diagonal  line.  It  has  N identical  singular  values.  But,  an 
array  with  a single  pixel  vertical  line  possesses  one 


10  1 0 1 


u i j 


10  10 
0 10  1 


10  10 
0 10  1 


1 0 
0 1 


1 0 
0 ] 


(9) 
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non-zero  singular  value. 


An 

SVD  is 
unitary 


interesting  and  useful  deterministic  property  of  the 
its  invariance  to  unitary  transforms.  Consider  a 
transformation 


— -c— R 


(10) 


Vlfrp 

where  LrL  = L_L  =1^.  It  is  easily  demonstrated  that  the 
characteristic  equations  of  11/  and  FF1  are  identical,  and 
hence,  the  singular  values  of  F and  Z are  the  same.  The 
eigenvector  matrices  are  related  by 


uz  - LcUf 


(11a) 


T 

-R-F 


(lib) 


As  a consequence  of  this  property,  it  is  observed  that  the 
singular  values  of  an  image  block  and  its  two-dimensional 
Fourier  or  Hadamard  transform,  for  example,  are  identical. 
Certainly,  the  deterministic  structure  of  such  arrays  is 
substantially  different. 


The  conclusion  of  these  observations  is  that  the 
singular  value  distribution  is  not  well  suited  for  the 
characterization  of  deterministic  line,  spot,  or  shape 
structure.  But  rather,  the  SVD  does  provide  an  indication 
of  the  structural  dependence  between  rows  and  columns  of  a 
pixel  array.  This  property,  as  will  be  shown,  is  extremely 
important  for  the  characterization  of  texture-like  regions 
of  an  image. 


Stochastic  Properties 

If  an  image  block  F is  considered  to  be  a sample  of  a 
two-dimensional  random  process,  then  the  singular  values 
s(n),  as  defined  by  eq. (1) , will  be  random  variables  since 
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the  SVD  is  a linear  transformation  of  F to  S.  Let  F be 

2 — — 

column  scanned  to  produce  the  N x 1 vector  f [4].  Then, 
the  image  block  may  be  stochastically  defined  by  its  joint 
probability  density  f (f ) . It  is  also  useful  to  describe  the 
stochastic  nature  of  the  image  block  in  terms  of  its  mean 
vector  and  its  covariance  matrix  K ^ . With  this 

information  available,  the  joint  density  of  the  singular 
values  P [S (1 ) ,S (2) , . . . ,S (N) ] and  the  singular  value  moments 
can  be  found,  in  theory,  however,  in  practice,  derivation  is 
d i f f icul t . 


Consider  first,  calculation  of  the  moments  of  S. 
Referring  to  eqs. (2)  and  (4),  the  mean  matrix  of  the  image 
block  can  be  written  as 

R 

E{F}  = M-,  = E{USVT)  = £ E{S(n)  u vj)  (12) 

- _F  n=l 


Similarly,  from 

eq. (3) , the 

means  of  the 

symmetric  matrix 

products  are 

T 

F.  { FF  ) 

= E{UAUT) 

(13a) 

e{ftf) 

= E{7Av'1  ) 

(13b) 

with  reference  to  eq.{6), 
be  expressed  as 

E{FFT}  = 


e(ftf}  = 


the  left-hand  terms  of  eq.(13)  can 


N 

N 

E^f  rT) 

n=l  ~n-n 

N 

xT 

II 

(14a) 

hT) 

rf=l  "n-n 

■ 

n = l 

(14b) 

where  R„  and  R_  denote  the  n-th  column  and  row  correlation 
-Cn  -Rn 

matrices  of  F,  respectively.  Continuing  one  step  further, 
the  means  of  the  symmetric  matrix  products  about  the  mean  of 
F are  given  by 
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(15a) 


E{ 


(F-M  )(F-M_,)T}  = 53 


n=l 

N 


E{(F-Mf)T(7-M?)}  = £*Rn 


(15b) 


n=l 


where  K and  K are  column  and  row  covariance  matrices  of 
Ln  Kn 

F. 


The  preceding  equations  provide  a link 


N 

N 

E{UAUT} 

- E 

«cn  = 

£ 

-Cn 

+ 

(16a) 

n=l 

n=l 

e{vavt) 

2W 

ii 

«Rn  = 

N 

E 

— Rn 

+ MpMp 

(16b) 

n=l 

n=  1 

between  the  mean  values  of  the  eigen-expansions  of  F and  its 
first  and  second  order  moments.  A special  case  of  practical 
interest  occurs  when  the  mean  matrix  M„  is  constant  valued 
and  the  row  and  column  covariance  matrices  are  identical 
from  row-to-row  and  column-to-column . For  this  case, 

E{UAUT)  = NK^  + MpMp  (17a) 

T T (17b) 
E{ VA v ) = NKr  + MpMp 

where  K ^ and  K R denote  the  column  and  row  covariance 
matrices,  respectively.  Thus,  the  mean  eigen-expansions 
assume  the  same  forms  as  the  array  column  and  row 
correlation  matrices.  Unfortunately,  it  is  not  possible  to 
proceed  further  with  this  approach  and  isolate  the  mean  of  A 
(equivalently,  the  second  moment  of  S)  from  eqs. (16)  or  (17) 
since  U and  V are  dependent  upon  A. 


The  moments  of  S can  be  found,  in  principle,  by  direct 
computation  of  the  means  of  the  eigenvalues  of  the 
characteristic  equation 
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c( A ) = det[FF  - X 


N 

FFT  - XI]  = Z C 

* r\  J- 


(18) 


1 = 0 


where  C . is  a function  of  the  elements  of  F. 

1 — 


roots  of  the  characteristic 
roots  of  eq.(17)  in  terms 


The  singular 

values  of  F are  the  positive  square  roots  of  the  eigenvalue 

equation.  Isolation  of  the 
of  F and  subsequent  moment 
calculation  has  proved  complex  even  for  general  2x2  pixel 
arrays.  Extension  of  the  approach  to  large  size  arrays 
seems  to  be  totally  intractable.  However,  there  is  some 

benefit  to  be  derived  from  consideration  of  the 

characteristic  equation  from  a stochastic  viewpoint. 

Consider  the  case  in  which  F is  a sample  from  an  array 

of  correlated  wide  sense  stationary  random  variables  with 

o 

common  variance  a and  a separable  covariance  matrix 
K = Kp  ® K , where  K and  K represent  column  and  row 
covariance  matrices,  respectively.  Now,  let  a 

Karhunen-Loeve  transform  be  performed  on  columns  and  rows  of 
F to  produce  the  random  array 


where 


Z = L^FLr 


^R  - — R 


(19) 

(20a) 

(20b) 


1 


and  L , L , £ , £ are  eigenvector  and  eigenvalue  matrices 
C R C R 


of 


and 


-c 

established 

invariant 


-R' 


as  noted.  From  eq. (11) , it  has  been 
that  the  singular  values  of  a matrix  are 
to  a two-dimensional  unitary  transformation 
because  the  characteristic  equation  is  not  changed  by  such  a 
transformation.  Since  the  K-L  transform  is  unitary,  it  is 
concluded  that  the  singular  values  of  the  correlated  array  F 
and  the  uncorrelated  array  Z are  identical.  This  may  seem 
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to  indicate  that  the  SVD  is  useless  for  characterizing  the 

structure  of  correlated  arrays.  But,  this  is  not  so.  The 

spatial  correlation  information  of  F has  simply  been 

converted  to  another  form  - the  variance  distribution  of 

To  show  this,  let  V „ be  an  N x N matrix  of  the  variance 

— a 

terms  of  Z in  spatial  correspondence  with  Z.  Then,  it  is 
easily  shown  that 

Y = I IT 

“Z  C R <21> 

where  the  vectors  y_  and  y contain  the  diagonal  terms  of 

L»  Iv  L 

and  respectively.  For  some  well-defined  processes,  such 

as  the  first  order  Markov  process,  can  be  determined 

analytically.  Thus,  to  summarize,  if  F is  wide  sense 
stationary  with  a separable  covariance  matrix,  then  the 
moments  of  F can  be  found  indirectly  from  the  moments  of  Z 
which  is  an  uncorrelated  array  with  an  energy  distribution 
given  by  eq. (21 ) . 

Attention  will  now  be  directed  toward  the  establishment 
of  the  joint  probability  density  of  the  singular  values  of 
F,  and  from  this,  the  singular  value  moments.  Consider  the 
outer  product  expansion 

N 

G = (F  - Mp)(F  - mf)t  = E (fn  - nn)(fn  - nn)T  (22) 

of  F about  its  mean  M „ where  f denotes  the  n-th  column  of  F 
- -F  -n  - 

and  n is  its  mean  (the  m-th  column  of  Mr) . Suppose  that  F 
n * 

is  considered  as  a collection  of  N independently  generated 

N-dimensional  vectors  f each  with  mean  n but  common 

-n  -n 

covariance  matrix  K^,.  Then,  G is  equivalent  to  the  well 
known  scatter  matrix  of  multivariate  statistical  analysis 
[5, p.547].  The  probability  density  of  the  scatter  matrix 
has  been  established  for  Gaussian  random  vectors.  In  this 
case 
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P(G)  = W(G,N,N,K_ -1) 

■“  — L 


(23) 


where 


WCG.K.N.K"1)  5 


N 1 

I K^1 1 2 I G I 2 


( N-K-l ) 


K K 


exp<-? 


^n(!)r(¥).  • • rM 


n.  N 


(24) 


is  called  the  Wishart  density  [6]  and  r(*)  is  the  gamma 
function. 


The  next  step  is  to  determine  the  joint  probability 

density  of  the  eigenvalues  of  G.  In  the  late  thirties 

several  mathematicians  found  this  density  for  a scatter 

2 

matrix  defined  about  its  experimental  mean  with  Kq  = o I 
[5, p.568].  From  these  results  the  eigenvalues  of  G are 
found  to  possess  a distribution  with  an  element  probability 


P(A)  dA ( 1 ) . . .dA(N)  = 


N 

n2 

' N 

_1  N i 

2 i>^=1  [ A ( .j  ) - A ( i ) ] exp  | 

. 2°2  i?ix(1)| 

N2  N 

(202)2  A 

[mj 

2 

dA( 1) . . .dA(N) 
(25) 


for  0 < A(n)  < ...  < A(l)  < 00 . The  second  moment  of  the 
n-th  singular  value  can  then  be  found  by  direct  evaluation 
of  the  integral 

Jr 00  /•  A ( 1 ) rA(N-l) 

[ ...I  X(n)  P(A)  dX(N)...dX(2)  dX(l)  (26) 

0 ^0  ^0 

No  simple  closed  form  solution  for  the  general  form  of  this 
integral  equation  has  been  found  so  far.  But  for  the  2x2 
case,  the  closed  form  solution  has  been  derived  as  the 
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following : 


E(S2 ( 1 ) } = a~ (2  + J)  (27) 
E(S2(2)}  = a2(2  - J) 

In  1960  James  [7]  generalized  the  solution  for  the 
distribution  of  eigenvalues  of  a scatter  matrix  for  an 
arbitrary  column  population  covariance  matrix  . The 
result  is  a complicated  functional  of  zonal  polynomials. 
Nevertheless,  the  distribution  can  be  numerically  evaluated, 
and  its  second  moments  can  be  determined  in  a similar 
fashion  to  eq. (26) . Then  by  the  d iagonal ization  procedure 
of  eqs. (19)  to  (21)  , it  is  possible  to  relate  the  second 
moment  of  the  singular  values  to  the  degree  of  spatial 
correlation  of  F for  a wide  sense  stationary  process  with  a 
separable  covariance  function. 

Figure  2 contains  a plot  of  the  singular  value  variance 
as  a function  of  the  spatial  correlation  factor  p for  a zero 
mean,  unit  variance  Gaussian  process  with  a separable 
covariance  matrix  K = K ® K_,  where  K_  = KD= 

These  theoretical  results  support  the  previous  contention 
that  the  singular  value  distribution  tends  toward  uniformity 
as  p approaches  zero,  and  toward  skewness  for  p large. 

Exper imental  Results 

Figure  3 contains  two  images  of  natural  texture,  grass 
and  ivy,  along  with  two  artificially  generated  fields 
obtained  by  convolutional  processing  of  two-dimensional 
fields  of  randomly  generated  pixels.  The  subjective  match 

between  the  natural  and  artificial  fields  appears  to  be  in 

✓ 

reasonable  agreement.  Figures  4 and  5 contain  plots  of  the 
singular  values  extracted  from  16  x 16  pixel  blocks  in  the 
center  of  each  image.  The  distributions  of  the  natural 
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(a)  natural  grass 


(b)  artificial  grass 


(c)  natural  ivy 


(d)  artificial  ivy 


Figure  3 . Examples  of  Natural  and  Artificial  Texture. 
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Figure  4.  Singular  Values  of  Natural  and  Artificial  Grass' 
Texture . 
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grass  and  ivy  are  seen  to  be  significantly  different,  but  the 
distributions  of  the  natural  and  artificial  pairs  are  quite 
close . 
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3 . Image  Processing  Projects 


The  image  processing  projects  described  in  the 
succeeding  pages  represent  a broad  variety  of  research 
topics  in  which  more  fundamental  signal  processing 
principles  are  utilized.  The  first  contribution  presents  a 
new  way  of  generating  binary  phase  holograms  via  the  digital 
computer.  Such  holograms  are  then  recorded  on  hardcopy 
devices  (see  Section  5)  and  after  film  development  are 
played  back  in  the  Institute  optical  laboratories.  The 
second  contribution  in  this  section  presents  rather 
encouraging  results  on  the  blind  deconvolution  process  of  a 
posteriori  restoration.  This  project  is  currently  coming  to 
fruition  with  quite  improved  results  for  arbitrary  space 
invariant  distortions.  A radar  imaging  task  is  next 
presented  in  which  multi-frequency  radar  returns  are 
coherently  processed  to  form  images  of  aircraft  mounted  on  a 
turntable.  Previous  reports  included  the  analysis 
associated  with  this  theory  and  this  report  presents  the 
images  resulting  from  such  processing.  This  section  closes 
with  a brief  summary  of  perceptual  model  color  image  coding. 
Results  are  presented  in  which  24:1,  48:1,  and  96:1 
compressions  are  obtained  with  good  visual  color  quality 
maintained . 


3.1  Double  Phase  Holograms,  A New  Way  of  Generating  Binary 
Holograms  (Supported  by  NSF  Grant  ENG-76-15318  and  AFOSR 
under  Contract  AFOSR-77-3285) 

Chung-Kai  Hsueh  and  Alexander  A.  Sawchuk 


Of  all  the  classes  of  computer  generated  holograms, 
binary  holograms  are  the  easiest  to  make.  Having  the 
Fourier  transform  of  the  object,  only  a plotter  or  a binary 
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display  device  and  photo-reduction  equipment  are  required  to 
make  a binary  hologram.  They  also  have  the  advantage  of 
superior  signal-to-noise  ratio  due  to  the  binary 
transmission  nature  of  the  hologram.  The  efficiency  of  a 
binary  hologram  can  also  be  increased  by  bleaching. 

The  major  difficulty  with  digital  holography  is  in 
finding  a method  to  record  complex  transmission  values  on  a 
photographic  or  other  optical  recording  medium  which  can 
store  non-negative  real  values.  Area  modulation  and  grey 
levels  are  used  to  represent  amplitude,  and  various  methods 
of  encoding  the  phase  by  position  modulation  [1]  within  a 
Fourier  transform  resolution  cell,  or  by  the  use  of  subcells 
[2-3]  have  been  developed. 


Since  any  complex  value  can  be  expressed  by  the 
summation  of  non-negative  multiples  of  +1,  -1,  +j  and  -j  Lee 
[2]  divided  each  cell  into  four  subcells  with  grey  levels  in 
each  subcell  to  represent  the  amplitude  of  each  vector 
component.  Burckhardt  [3]  further  found  that  instead  of 
using  four  vectors  to  represent  a complex  value,  three 
vectors  pointing  in  the  directions  0,  2tt/3  and  -2ir/3  can  be 
used  to  represent  any  complex  value. 


The  Double  Phase  Hologram  (DPH)  introduced  here  is 

based  on  the  principle  that  any  complex  transmission  value 
1 6 

AeJ  inside  the  unit  circle  (0  < A < 1)  can  be  decomposed 
into  the  sum  of  two  constant  magnitude  vectors  as  expressed 
by  [4] 


H = Ae 


je  = 11(0+*)  + (0-^) 


where 


= cos  A, 


0 < <Ji  < tt/2 


-91- 


In  this  expression  we  have  dropped  the  (u,v)  coordinates  of 
the  transform  domain  for  simplicity.  The  corresponding 
coordinates  in  the  image  plane  will  be  designated  (x,y) . 
The  two  terms  on  the  right  side  of  Eq.  (1)  are  designated  H 
and  H respectively,  and  these  are  shown  in  Fig.  1.  The 
decomposition  can  be  applied  to  represent  a complex  value  in 
a binary  hologram  by  using  the  two  phase  quantities  instead 
of  the  three  non-negative  quantities  in  a Burckhardt 
hologram. 

These  two  phase  quantities  are  implemented  by  detour 
phase  as  shown  in  Fig.  2.  Each  Fourier  component  on  the 
hologram  consists  of  two  subcells  with  each  subcell  having 
an  open  slit  of  fixed  area  in  the  opaque  background.  The 
position  is  determined  by  the  decomposed  phases  from 

dl  = d(e+i|0/2Tr  (3) 

and 

d2  = d(e-ip)/ 2ir  (4) 

where  d is  the  size  of  one  Fourier  resolution  cell.  If  we 
consider  the  subcells  as  a whole,  then  each  component  has 
the  correct  complex  value  in  the  direction  of  the  rirst 
diffraction  order  along  the  v direction.  Therafc^e  ‘ ’c 
image  appears  at  the  first  diffraction  ordt'r  in  th?  na;  k 
focal  plane  of  the  lens.  This  type  of  hologram  3 ^<*1’  id  a 

Double  Phase  Hologram  (DPH) . 


From  Eq.  (2)  the  angle  \p  is  related  to  the  amplitude  A 


by 


A = cosip 


0 < ip  < 7t/2 


(5) 
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and  its  derivative  is  given  by 


r 


= sinij/  0 < < tt/2  (6) 

A 

For  the  DPH  is  quantized  to  a certain  level  iji  and  the 
amplitude  becomes 

A A 

A = cosi|>  (7) 

rather  than  A.  As  shown  in  Eq.  (6),  the  derivative  of  A is 
an  increasing  function  in  the  defined  range.  Therefore  the 

A 

same  amount  of  quantization  error  in  ip  results  in  smaller 

/A 

error  in  A when  A is  close  to  1 (ip  close  to  0)  and  larger 

A 

error  in  A when  A is  close  to  0 (ip  close  to  it/2)  . This 

indicates  that  it  is  better  to  have  an  uniform  spectrum  in 

order  to  reduce  the  quantization  error.  An  object  with  zero 

phase  usually  has  a very  large  dynamic  range  and  energy  is 

concentrated  at  low  frequency  components.  To  spread  out  the 

spectrum  we  can  either  impose  a random  phase  on  the  object 
or  use  deterministic  diffusers  [5].  Other  iterative  phase 
coding  techniques  [6-7]  can  also  be  used  to  smooth  out  the 
spectrum  in  a desired  form  and  therefore  reduce  the 

quantization  error  in  the  DPH. 

It  has  been  found  that  although  the  phase  coding 
techniques  reduce  the  dynamic  range  of  the  Fourier 

transform,  they  also  introduce  speckles  due  to  aliasing 

error.  The  space-bandwidth  iterative  method  with 

constrained  bandwidth  [8]  has  been  used  to  reduce  this 

noise.  Here  we  can  use  a similar  idea  to  implement  a DPH. 
First  of  all  a Fourier  transform  with  constrained  bandwidth 
is  obtained  by  the  iterative  method.  Then  only  the  center 
part  of  the  transform  is  used  in  plotting  a hologram.  By 
doing  this  we  can  reduce  the  speckle  noise  while  preserving 
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the  same  space-bandwidth  product.  The  uniform  transform 
required  for  correct  phase  quantization  in  a DPH  is  also 
maintained . 


. 

L 


The  implementation  of  a DPH  is  based  on  the 
decomposition  of  a vector  into  two  other  constant  vectors  as 
given  by  Eqs.  (1)  and  (2) . By  defining 


H + P = 2HX 


(8) 


and 


H - P = 2H2 


(9) 


where  P is  a vector  called  the  parity  term  shown  in  Fig.  1, 
the  decomposition  of  Eqs.  (1)  and  (2)  can  be  rewritten  as 
[9] 


H = Hj+H2  = £(H+P)  + J(H-P) 


(10) 


If  and  H2  were  overlapped  at  the  same  position,  P and  -P 
will  cancel  out  and  the  correct  transform  is  obtained. 
However  due  to  the  method  of  implementing  the  DPH,  these  two 
terms  do  not  exactly  cancel  out,  generally  resulting  in 
parity  noise. 

Let  us  denote  h(x,y)  and  p(x,y)  as  the  inverse  Fourier 
transforms  of  H(u,v)  and  P(u,v)  respectively.  Using  the 
configuration  in  Fig.  2,  it  can  be  shown  that  the 
reconstruction  at  the  first  diffraction  order  along  the  y 
direction  is  given  by 

0(x,y)  = hs(x.y)  + PsU.y)  (11) 


with 
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J 


(12) 


h (x,y)  = d Wh(x,y)  [h(x,y)*coTnb(dx,  dy)  ] 


o 

ps(x,y)  = jd  Wp (x,y) [p (x,y)*comb (dx, dy) ] 


(13) 


where 


Wh(x,y)  = sine  (dx,  dy) 


(14) 


ard 


W (x,y)  = -sinc(dx/2  , dy)sin(iTdx/2) 

P (lb) 

are  the  weighting  functions  for  the  desired  function  and  the 
noise  term. 

Both  W, (x,y)  and  W (x,y)  are  drawn  in  Figs.  3 and  4 
h p 

respectively.  At  x = 0 the  noise  term  p(x,y)  is  strongly 
attenuated  by  sin(irdx/2).  Since  h(x,y)  is  space-limited  we 
may  assume  that  p(x,y)  is  also  space-limited  to  a small 
region.  Therefore  in  the  region  of  interest  (zero  order) 
noise  is  less  noticeable.  On  the  other  hand  at  x = 1/d 
(first  order)  sin(Ttdx/2)  reaches  its  maximum  and  sinc(dx/2) 
does  not  fall  off  too  much  while  h(x,y)  is  strongly 
attenuated  by  sinc(dx)  at  x = 1/d.  Therefore  the  error  term 
dominates  at  the  first  diffraction  order. 

Each  cell  can  be  further  divided  into  more  subcells  as 
long  as  they  sum  up  with  the  correct  transform  value. 
Suppose  there  are  2n  subcells  in  each  cell  with  and  H2 

plotted  alternatively.  We  can  proceed  as  before  and  obtain 
Eqs.  (11)— (14)  with 


-96- 


Wp(x,y)  = — pp— sinc(dx/2n,  dy)  {sin<J>-sin3<)>+.  . , + (-l)n  ^sin(2n-l)<|>}  (16) 

where  <J>  = ■rrdx^n.  The  curves  of  Wp  (x,y)  along  the  x 
direction  for  n = 2,  4,  8 are  shown  in  Figs.  5-7 
respectively. 

Some  conclusions  can  be  made  based  on  the  above 
observations : 

(1) .  As  n increases  the  noise  term  around  x - 0 is 
expected  to  reduce  and  we  can  get  better  reconstruction  by 
using  more  subcells. 

(2)  . The  dominating  noise  terms  appear  at  the  +n 
orders . 

(3) .  Other  noise  term  orders  are  strongly  attenuated 
and  most  of  them  are  on  the  zeroes  of  the  weighting 
function . 

We  have  performed  experiments  to  verify  the  above 
analysis.  The  letter  'P'  was  used  as  the  object  and  was 
precompensated  according  to  the  weighting  function  W (x,y) 
shown  in  Fig.  3.  Fig.  8 shows  the  reconstruction  with 
random  phase.  Speckle  due  to  aliasing  error  is  clearly 
observed.  Fig.  9 shows  the  reconstruction  when  the 
space-transform  iterative  method  with  constrained  bandwidth 
is  used  and  only  the  central  part  is  plotted.  The  image 
quality  is  considerably  improved. 
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3.2  A Technique  ot  A Posteriori  Restoration — Results  of  a 
Computer  Simulation 


John  Morton 


Illustrated  in 
point-spread-f unction 
reported  herein. 


Figure  1 is  the  blurring 

(PSF)  used  in  the  computer  simulation 


Figure  2 compares  the  log^Q  of  the  magnitude  of  the 
optical  transfer  function  (OTF)  with  the  corresponding 
estimate.  For  comparison  purposes  values  of  the  log  of 
the  magnitude  below  2 have  been  set  to  2.  The  estimate  used 
the  method  of  Cannon  [1]. 


Figure  3 compares  the  phase  of  the  OTF  with  two 
estimates.  Both  estimates  used  the  method  as  previously 
reported  [2-4];  the  differences  in  the  two  estimates  is  that 
Estimate  1 relaxed  one  assumption  implicit  in  the  method, 
whereby  Estimate  2 did  not.  The  comparisons  for  the  phase 
of  the  OTF  on  the  frequency  axes  are  given  in  Figure  4. 

Figure  5 presents  the  image  degraded  by  the  PSF  in 
Figure  1 together  with  restorations  assuming  different  a 
priori  knowledge.  The  assumptions  underlying  the 
restorations  are  given  in  Table  1. 

Upon  comparison  of  the  magnitude  of  the  optical 
transfer  function  and  the  estimate  of  the  magnitude,  it  is 
clear  that  the  estimate  adequately  retains  the  essential 
features . 

Upon  comparison  of  the  phase  of  the  optical  transfer 
function  and  the  estimates  of  the  phase,  it  is  evident  that 
phase  Estimate  1 is  reasonable  near  the  axes  and  not  as 
accurate  at  a distance  from  the  axes.  Because  images  in 
general  have  little  energy  concentration  in  those  frequency 
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(d)  magnitude  of  OTF  given 
phase  estimate  = Estim 


(c)  magnitude  of  OTF  given 

phase  of  OTF  estimated  as  0 


(e)  magnitude  of  OTF  estimated  (f)  magnitude  of  OTF  estimated 

phase  of  OTF  given  phase  of  OTF  estimated  as  0 


(g)  magnitude  of  OTF  estimated 
phase  estimate  = Estimate  2 


(h)  magnitude  of  OTF  estimated 
phase  estimate  = Estimate  2 


Figure  5 (continued) 
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Figure 

OTF 

OTF 

Power  Spectrum 
of 

magnitude 

phase 

Undegraded  image 

b 

given 

given 

given 

c 

given 

estimated  as  0 

given 

d 

given 

Estimate  2 

given 

e 

estimated 

given 

estimated 

f 

estimated 

estimated  as  0 

estimated 

g 

estimated 

Estimate  2 

estimated 

h 

estimated 

Estimate  1 

estimated 

Table  1.  Key  to  restorations. 
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components  away  from  the  axes,  the  poorer  results  away  from 
the  axes  are  not  very  important  from  a restoration  point  of 
view. 

From  the  phase  comparisons  it  may  be  noted  that  phase 
Estimate  2 is  not  particularly  good. 

Comparing  restorations,  it  is  apparent  that  the 
estimate  of  the  magnitude  of  the  OTF  is  reasonable  to  the 
extent  of  achieving  an  adequate  restoration.  In  fact,  the 
visual  quality  of  the  restoration  using  the  estimate  of  the 
magnitude  of  the  OTF  seems  to  be  superior  to  the  restoration 
given  knowledge  of  the  magnitude  of  the  OTF.  In  general, 
upon  restoration,  the  frequencies  in  the  vicinity  of  the 
zeros  of  the  OTF  will  be  magnified  the  most  by  the  Wiener 
filter.  If  one  knows  the  magnitude  of  the  OTF  exactly,  at 
times  the  restoration  will  contain  periodic  artifacts 
corresponding  to  those  frequencies  that  were  boosted  the 
most  by  the  filter.  On  the  other  hand,  the  inaccuracies 
inherent  in  the  estimate  of  the  magnitude  of  the  OTF  cause 
the  frequencies  in  the  vicinity  of  the  zeros  to  be  boosted 
less  than  the  boost  assuming  knowledge  of  the  magnitude  of 
the  OTF.  As  a result,  when  assuming  the  estimate  of  the 
magnitude  of  the  OTF,  the  periodic  artifacts  are  missing. 

Although  the  visual  quality  of  the  restoration  is 
usually  superior  when  using  the  estimate  of  the  magnitude  of 
the  OTF,  in  terms  of  minimum  mean  squared  error,  restoration 
b is  superior. 

Upon  comparing  restorations  c and  d and  restorations  f 
and  g,  it  is  generally  agreed  that  Estimate  2 of  the  phase 
is  not  significantly  better  with  respect  to  restoring  than 
an  estimate  of  zero  for  the  phase. 
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The  superiority  of  Estimate  1 of  the  phase  is  evident. 
For  example,  note  less  distortion  in  the  face  of  the  man  in 
h opposed  to  f and  g.  Also  note  the  streaks  of  light  in  the 
pants  of  the  man  and  across  the  chest  and  below  the  knees  of 
the  woman  in  f and  g.  These  anomalies  are  absent  in  h. 
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3.3  Turntable  Radar  Imaging 

Chung  Ching  Chen  and  Harry  C.  Andrews 

Introduction 

Radar  systems  use  active  device  to  radiate  radio 
frequency  signals  and  receive  their  returns.  Various 
information  about  the  targets  can  be  inferred  by  processing 
those  echoed  data  [1], 
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In  two  dimensional  radar  imaging  system  the  two 
geometric  coordinates  are  usually  range  and  azimuth  [2]. 
Range  is  the  direction  along  which  the  signal  is  transmitted 
and  received.  Azimuth  is  the  direction  orthogonal  to  the 
range  direction  and  along  the  relative  notion  of  radar  and 
targets.  The  range  information  is  obtained  by  timing  the 
radar  return  by  either  short  pulse  or  pulse  compression 
techniques.  In  the  clasrical  radar  systems  azimuth 
information  is  provided  by  using  very  long  antenna  whose 
illuminating  pattern  is  very  narrow  along  azimuth,  because 

0 = X/L  (1) 

where  6 is  the  effective  angle  of  antenna  pattern,  A is  the 
wavelength  and  L is  the  length  of  the  antenna.  To  achieve 
small  0 comparable  to  the  range  resolution  usually  requires 
impractically  long  antennas.  Coherent  radar  systems  get 
around  this  problem  by  synthesizing  without  actually 
implementing  long  antennas  created  by  the  motion  of  the 
targets  and  the  radar  [3],  [4],  [5], 

In  this  paper  we  study  a radar  system  which  has 
mathematical  similarity  with  the  computer  aided  tomography 
(CAT)  used  in  medical  applications.  Degrees  of  freedom 
(DOF)  (6]  of  such  a system  is  derived  for  the  purpose  of 
efficient  computation  and  minimal  data  storage  requirement. 
Physical  differences  between  the  two  systems  are  analyzed 
and  incorporated  in  the  tailored  reconstruction  algorithm 
unique  to  the  underlying  radar  system  to  reconstruct  the 
target  images. 

Experiments  are  done  in  support  of  the  theory  developed 
in  this  work. 
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The  RAT  SCAT  Facility  - Data  Acquisition 

The  RAT  SCAT  (RAdar  Target  SCATter  Site)  [7]  facility 
is  used  to  obtain  the  data  for  our  processing.  A target  to 
be  imaged  is  placed  on  a rotator  at  a distance  r^  from  the 
radar  to  the  rotation  center  as  in  Fig.  1.  The  reference 
sphere  S sits  at  distances  r^  from  the  radar  R and  r^  from 
the  rotation  center  C.  As  shown  in  Fig.  2 there  are  two 
rectangular  coordinates  systems.  (£»n)  and  (x,y)  , with  the 
later  system  making  an  angle  0 from  the  former  one.  We 
assume  that  (£,n ) is  fixed  with  the  target  and  (x,y)  is 
fixed  with  the  radar  or  ground.  The  radar  radiates  the  same 
set  of  step  frequencies  { f^}  at  discrete  aspect  angles  from 
the  set  {6^}.  Normally  the  step  frequency  Af  and  step 
single  A0  are  fixed,  and  we  assume  so  in  our  case.  The 
local  oscillator  at  S takes  the  signal  from  R to  S and  beats 
it  with  the  signal  reflected  from  the  target  to  form 
in-phase  and  quadratic-phase  data  components. 

Target  Reflectivity  Function 

Define  the  reflectivity  function  f (£  ,n ) of  the  target 
as  the  ratio  of  the  received  signal  due  to  the  target  point 
at  (C  , n)  and  the  transmitted  signal,  as  in  Fig.  2.  The 
metallic  target  looks  specular  to  the  radar  because  the 
wavelength  A used  is  very  small  compared  to  the  target  size. 
There  is  also  shadowing  effect  because  of  the  non  convexity 
of  the  target  surface.  Thus  f (£ ,n ) is  actually  a function 
of  0.  However  we'll  assume  for  the  moment  that  f ({;  ,n ) is 
independent  of  0 and  we'll  derive  the  PSF  as  well  as  DOF  of 
such  a system  which  offer  much  insight  into  the  mathematical 
properties  of  the  systems. 

PSF  Of  The  Imaging  System 

In  Fig.  1 let's  assume  that  angles  a and  g are 
virtually  zero.  The  distance  between  the  target  point  (£ ,n ) 
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on  rotator/ 

h 

ro 

- 

Fig.  1. 

relation  among  radar, 
reference  sphere 

target  and 

and  the  radar  is 


r(r0,e,£,n)  - [(rQ-x) 2 + y2]% 

~ (rQ-x), 

where  we  have  assumed  r^x  >>  y for  all  target  points 
and 


(2) 

(x,y) 


x = £cos0  + nsin0 


(3) 


' y = ncos©  - £sin0 

Differentiating  (1)  with  respect  to  e and  incorporating 
Eq.  (2)  one  gets 


3r(ro,0,5,n) 

— = £sin0  - ncos0  = -y  (4) 

30 

where  r^  -x  % r(r  ,e,£,n)  has  been  assumed.  It  is  observed 
from  Eqs.  (2)  and  (4)  that  the  lines  of  constant  range  and 
the  lines  of  constant  doppler  are  parallel  to  the  y,x,  axis, 
respectively. 

Define 

a 

g(x,0)  = f f(e,n)dy 


where  a is  the  maximal  radial  extent  of  the  target,  and 

fk(t)  - A cos(2irfkt+<t>)  k = 1,2 M (6) 

which  are  the  signals  transmitted  from  the  radar.  The 
signals  returned  to  S from  the  target  line  mass  at  range 
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(7) 


z(x,i,k,t)  = BgCx.0^  cos^2nfk^t — 9_J 

where  the  signal  power,  propagation  decay  and  the 
reflectivity  phase  have  been  absorbed  into  the  complex 
constant  B.  The  total  return  is 


z(i,k, t) 


z(x,i,k,t)dx 

‘ -a 


= B Ja  g(x , 9.)  cos  ^2TTfk  (t- -2X)-*-»]  dx 


(8) 


The  signal  along  path  RS  received  by  S assumes  the  form 


r (k, t)  = C 


cos 


[2'fk(t-T:)+*] 


(9) 


Beating  signals  (7)  and  (8)  and  appropriately  filtering 
the  result  in  in-phase  and  quadratic-phase  components 

I(i,k)  - D | g(x,9i)cos  1 2 7T f k^^  + 27rfk  1 'c — “ * dx 


Q(i.k)  = -D  |a  g(x,9.)sinj2TTfk^  + Zirf^1-^0  ^ | dx 


(10) 


were  not  for  the  linear  phase  factor  Eq.(10)  would  be  the 

Fourier  compo^nts  of  the  shadow  gram  g(x,e^)  at  angle 

and  frequency  k,  (8),  [9],  [11]  and  would  take  values  on 

c 

the  ring  area  as  in  Fig.  3 just  like  that  in  the  tomography. 
DOF  of  the  tomographic  system  has  been  analyzed  in  detail  by 
McCaughey  and  Andrews  [10]  with  the  same  geometry  as  shown 
in  Fig.  4.  Because  of  the  mathematical  similarity  we'll 
apply  some  of  the  results  in  [10]  to  find  the  DOF  of  the  raw 
radar  data,  keeping  in  mind  the  physical  differences  such  as 
the  angle-dependence  of  reflectivity  function  mentioned 
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above . 


DOF  Of  The  Data 

The  DOF  analysis  tries  to  determine  the  number  of  truly 
independent  data  under  certain  system  geometry.  A 
conventional  way  to  find  the  DOF  is  to  find  the  correlation 
matrix  [r ] of  the  point  spread  function  of  the  system  and 
diagonize  it.  Roughly  speaking,  the  rank  of  [ r ] represents 
the  DOF. 


In  order  to  easily  apply  the  results  in  [10]  directly 
to  our  radar  data,  let's  normalize  the  spatial  extent  of  the 
target  as  in  [10].  If  the  physical  maximum  and  minimum 
radio  frequency  are  f[n^n  and  f^^  corresponding  to  the  2 way 
wavelengths 


■>  1 „ C _ 1. 5x10  m/sec 

max  " T~  x ? T~. 

min  mm 


— — * 1.5xl0^m/sec 


then  the  normalized  frequencies  will  occupy  the  spectral 
ring  from  _ a 


fi = r 


f = — 

2 X • 
min 


as  in  Fig.  5. 


From  [10]  the  (i,m,k,l)th  entry  r^’™  of  [T]  is 
rk’l  = | | exp[-j21Tuk(Cc°sei+nsin0i)] 

exo[+j  2ttu^  (CcosOjj+nsin^)  ]d£dn 

[211(0^-21^^003  (0i-0m)+u][) 

[uk'2ukuicos(9i-V+ui:|,i 
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Fig.  5.  spatial  domain  and  frequency  domain  of 
a spatially  normalized  target 


m are 


where  (i,k)  and  (m,l)  are  any  two  data  points;  i and 

azimuth  angle  indices  and  k and  1 are  frequency  indices 

along  radial  direction.  (•  ) is  the  first  order  Bessel 

function  of  the  first  kind,  R is  the  unit  circle  and  0 < 0., 

— l 

0m  < 2it  . f^  < 1^2*  0ne  can  rewrite  Eq.  (20)  as 

i,m  J]_[2rrp (k,  1 , i ,tn)  ] 
rk , 1 p (k , 1 , i , m)  ( 1 3 ) 

where  p(k , 1 , i ,m)  = [U^-2U^U^C°S  (0^-0^)  is  the  distance 

between  the  points  (i,k)  and  (m,l)  in  the  frequency  domain. 

Thus  one  can  think  of  the  data  as  samples  of  a stationary 

J-.  (2irp) 

field  with  correlation  function  , which  has  a power 

p 

spectrum  (Fourier  transform  of  the  correlation  function) 
circ(r)  with  radial  cutoff  frequency  1,  as  shown  in  Fig.  6. 
According  to  the  sampling  theorem,  in  order  to  avoid 
aliasing  in  the  spatial  domain  the  sampling  interval  in  the 
frequency  domain  must  be  < 1/2.  This  agrees  with  the 
intuitive  fact  that  to  adequately  represent  an  object  of 
maximum  spatial  extent  S,  the  sampling  interval  in  its 
frequency  domain  should  not  be  larger  than  1/S  (S=2  for  our 
normalized  target) . 

Actual  Reconstruction  Method 

Although  the  radar  system  has  a PSf  similar  to  that  of 
the  tomographic  system,  there  are  physical  differences  as 
mentioned  before  which  discourage  one  from  applying  the  same 
algorithms  as  used  in  the  tomographic  system  to  reconstruct 
the  target  images. 

In  Fig.  7 let  (x^,y^)  and  ,y^)  be  the  coordinates  of 
the  same  targets  at  two  aspect  angles  whose  difference  is  <j> . 
Similar  to  Eq. (3)  we  have 
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x2  = xiCOS(f)  + y^sin<j» 
y2  = yicos,«'  " x^sin4) 


(14) 


Thus  from  (x^y^)  to  (x2,y2)  the  range  of  the  point 
undergoes  a change  of 


Ar(xi>y1,x2>y2)  = x2~xl 

= x-^cos(J)+y^sin4i-x^ 

= x^(cos<|>  - l)4-y  ^sin4> 


(15) 


which  creates  different  phase  histories  to  different  points. 
If  A is  so  small  that 


cos<{>  ~ 1 


and 


Ar(x1,y1,x2,y2)  ss  y^in^  % yx4>  (16) 


which  is  independent  of  x^  or  the  range,  then  the  range 


processing  is  separable  from  the  azimuth  processing,  the 


doppler  frequency  is  which  is  a function  of  azimuth 


only.  Thus  the  range  and  azimuth  information  can  be 
resolved  by  taking  a 2-D  Fourier  transform  which  is 
separable  (13]  on  the  data  spanning  small  azimuthal  angle. 
To  insure  the  separability  the  range  walking  has  to  be 
prohibited  by  limiting  Ar  in  (15)  to  be  less  than  one  range 
bin  width. 


Fig.  8 shows  data  section  of  such  a small  fan  angle 


Let  B 

r 

and  B, 


B 


be  the  range  bandwidth  with  step  frequency  A f ^ and 
be  the  minimum  and  maximum  azimuth  bandwidths  at 


^ and  f2  respectively  with  step  frequencies  Afz 


f 

then 


and  A f „ 


I 

1 
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r 


Br  = M x A fr  = f2" f i 


BZl  * Na  x 4£Zl  * Na  x £140 


B 


* No  x Af  = N X f_A0 


a 


a 


(17) 


where  A0  is  the  azimuth  step  angle  in  radian,  M is  the 
number  of  samples  in  radial  frequencies  and  Ng  is  the  number 
of  azimuth  points  coherently  processed. 


If  R is  the  radial  spatial  extent  and  Z^ , Z2  are  the 
azimuth  spatial  extent  after  applying  the  discrete  Fourier 
transforms,  then 


R = 


The  dual  relations  are 


Ar  = W— 

Br 


Z,  = 


XT 


Zo  = 


1 


A F 


(18) 


Azl  = F 


Az2  " F 


(19) 


where  Ar,  AZ^,  and  AZ£  are  resolutions. 

Experimental  Results 

Utilizing  the  above  principles  actual  radar  returns 
were  processed.  A coherent  angle  of  6.4°  (equalling  32 
pulses  in  azimuth)  was  assumed  and  2-D  FFT  is  taken  over 
these  pulses.  The  results  are  presented  in  Fig.  9 for 
various  angle  of  rotation. 

To  improve  the  image  quality  nocoherent  integration  is 
performed  with  the  range  cross-range  images  as  in  figure  9. 
With  only  7 looks  noncoherently  summed  (at  30°  angle 
intervals)  the  image  of  figure  10a  results.  This  is  a 
considerable  improvement  and  clearly  shows  the  outline  of 
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a)  4 looks  0°-25 . 6° 
(nose) 


b)  4 looks  76.8°- 102. 4C 
(broadside) 


c)  4 looks  128°- 153. 6°  d)  4 looks  153 . 6°- 179 . 2° 

(tail) 


(The  Radar  is  positioned  on  the  right) 

Figure  9.  6.4°  coherence  in  Azimuth  at  various 

positions  of  rotation. 


a) 


7 looks  (6.4°  coherence) 
spaced  30°  apart 


b)  28  looks  (6.4°  coherence) 
spaced  6.4°  apart 


c)  56  looks  (3.2°  coherence) 
spaced  3.2°  apart 


d)  14  looks  (12.8°  coherence) 
spaced  12.8  apart 


j 


Figure  10.  F102  Imaged  for  various  parameters 
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the  characteristic  delta  wing  of  the  F102  aircraft.  By 
noncoherently  integrating  28  looks  one  obtains  the  results 
of  figure  10b  in  which  a more  clear  image  results.  To 
investigate  the  degree  of  coherence  necessary  (and  allowable 
before  "range  walking"  occurs)  figures  (10c)  and  (lOd) 
present  result  for  3.2°  coherence  angles  and  12.8°  coherence 
angles.  In  both  cases  the  aircraft  is  still  clearly  visible 
although  a certain  amount  of  degradation  is  beginning  to  be 
apparent  in  both  cases. 

A second  aircraft  was  imaged  using  the  same  parameters 
as  developed  above.  This  aircraft  was  an  F5E  and  is  shorter 
with  stubby  wings  and  wingtip  pontoons.  Figure  11  presents 
the  results  of  imaging  this  airframe  with  28  looks  at  a 
cross  range  coherence  integration  of  6.4°.  In  this  figure 
the  sensitivity  to  knowledge  of  the  center  of  motion 
(rotation  in  our  case)  is  investigated.  This  amounts  to 
adjusting  the  linear  phase  factors  described  earlier  . 
By  visual  inspection  probably  the  best  "focused"  image  is 
that  of  figure  lib.  The  advantage  of  using  the  techniques 
outlined  in  this  figure  is  that  iterative  focus  can  be 
developed  for  visual  inspection  on  the  stored  one 
dimensional  radar  returns. 

The  final  figure  (figure  12)  presents  a summary  of 
photographs  for  the  F5E  and  F102  airframes  for  both  azimuth 
and  elevation  plots.  Because  all  parameters  are  fixed  for 
these  images,  scales  are  preserved.  Consequently  it  is 
clear  that  the  F5E  is  a smaller  aircraft  and  naturally  has  a 
different  azimuth  and  elevation  projection  than  does  the 
F102. 
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a)  28  looks  (6.4°  coherence) 
centered  at  range  bin  80 


b)  28  looks  (6.4°  coherence) 
centered  at  range  bin  79.75 


c)  28  looks  (6.4°  coherence) 
centered  at  range  bin  79.5 


d)  28  looks  (6.4°  coherence) 
centered  at  range  bin  79 


Figure  11.  F5E  and  Center  of  Rotation 
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a)  F102  Azimuth  Image  b)  F102  Elevation  Image 

(28  looks)  (56  looks)  (1/2  scale) 


Conclusion 


This  paper  has  attempted  to  present  the  theory  of  high 
range  resolution  radar  imaging  from  both  a radar  systems 
viewpoint  and  a degrees  of  freedom  or  numerical  analysis 
viewpoint.  Similarity  with  the  computer  aided  tomographic 
scanner  imaging  technology  is  pointed  out.  However  the 
differences  between  the  two  systems  are  emphasized  and  a 
radar  unique  reconstruction  algorithm  is  developed  for 
combined  coherent  and  noncoherent  imaging.  The  actual 
reconstruction  method  is  explained  and  experimental  results 
developed  to  illustrate  the  theories  presented.  The 
pictorial  images  resulting  from  the  computational  procedures 
are  surprisingly  recognizable  and  suggest  that  these 
techniques  may  have  some  practical  application  in  the 
future . 
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3.4  Perceptual  Model  Coding  (Supported  by  WPAFB  under 
Contract  F--33615-77-C-1016 ) 

Charles  F.  Hall  and  Harry  C.  Andrews 


The  results  of  the  human  visual  system  perceptual  model 
for  color  image  compression  represent  coding  results  which 
are  indeed  encouraging.  The  model  takes  into  account  the 
nonlinearities  of  the  brightness  response  of  the  eye  as  well 
as  the  necessary  spectral  and  spatial  processes  to  describe 
non-temporal  visual  phenomena  in  a mathematical  formulation. 
Figure  1 presents  the  basic  formulation  utilized  in  this 
work.  Tristimulus  values  of  an  image  (R  (x , y)  ,G  (x ,y)  ,B  (x, y)  ) 
are  passed  through  a linear  transformation  T for 
luminous-like  and  chrominance-like  signals  (d^,  d£,  and  d^ 
respectively) . The  nonlinear  logarithm  device  provides 
large  dynamic  range  sensitivity  and  the  summing  functions 
thereafter  provide  color  constancy  on  the  chrominance 
channels  e2  and  e^.  The  gain  constants  C2»  develop  a 
color  space  where  just  noticeable  color  differences  (jncd's) 
are  unit  circles.  Finally  the  spatial  filters  h^(x,y), 
h2 (x,y) , and  h3(x,y)  remove  those  frequency  aspects  of  an 
image  which  the  eye  does  not  see  as  well  as  emphasizing 
those  aspects  to  which  the  eye  is  particularly  susceptible. 

This  model  has  been  used  with  considerable  success  in 
compressing  color  images  to  as  low  as  1,  1/2,  and  1/4  bit 
per  pixel.  The  color  plates  accompanying  this  section 
present  these  results  in  pictorial  form. 
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Figure  1.  Perceptual  Model  of  Color  Vision. 


Figure  7.9  Perceptual  Power  Spectrum  Coded  LAKE  Image  (N 


rum 


Figure  7.  11  Perceptual  Power  Spectrum  Coded  BUILDING  Image 
(N  = 512) 


\ V 


/ 


Figure  7.  12  Perceptual  Power  Spectrum  Coded  BABOON  Image 
(N  = 512) 


4.  Smart  Sensor  Projects 


This  section  reports  on  the  successes  of  the  smart 
sensor  CCD  chip  development  under  contract  at  the  Hughes 
Research  Laboratories.  Earlier  reports  presented 
descriptions  of  two  circuits  under  development,  the  Sobel 
circuit  and  Circuit  II  which  included  considerably  more 
complex  processes.  Both  circuits  have  been  demonstrated  to 
produce  proper  processing  results.  These  past  six  months 
have  seen  the  implementation  of  Circuit  II  in  real  time  TV 
(2  MHz  reduced  horizontal  bandwidth)  with  successful 
results.  It  is  anticipated  that  minor  redesign  will  achieve 
the  nominal  10  MHz  real  time  TV  bandwidth.  In  addition 
progress  is  being  made  on  other  fronts.  New  circuit  designs 
are  in  process  and  a two-year  plan  is  presented  in  the  last 
contribution  to  this  section. 


4.1  Charge  Coupled  Device  Technology  for  Smart  Sensors 


Graham  R.  Nudd,  Paul  A.  Nygaard,  and  Gary  D.  Thurmond 


Abstract 

This  paper  describes  our  continuing  work  II]  to  design, 
fabricate  and  test  charge  coupled  device  (CCD)  circuits  for 
image  preprocessing.  Two  test  chips  containing  six 
processing  algorithms  have  now  been  fabricated  and  tested. 
The  processing  functions  are  described  together  with  the 
circuit  implementation  and  a performance  evaluation. 

Processor  Development 

We  have  completed  the  design  and  fabrication  of  two 
test  chips  as  shown  in  Figures  la  and  lb.  These  circuits 
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Figure  la.  Photomicrograph  of  CCD  Sobel  Circuit. 
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are  two  phase  surface  channel  devices  with  8 micron  gate 
lengths.  N-type  silicon  is  used  to  achieve  maximum  speed. 
The  algorithms  implemented  are: 


Sobel  edge  detection: 

fs  = 1/8  { | (a+2b+c)  - (g+2h+i) | 
+ | (a+2d+g)  - (c+2f+i)|} 


(la) 


|sxl  + |Sy| 


(lb) 


Local  averaging: 


m 


1/9  (a+b+c+d+e+f+g+h+i) 


Unsharp  masking: 


usm 


(1-a)  e+af 


(2) 


(3) 


Binar ization: 


s-i: 


f_  < e 
ci  - 

f_>  e 
m 


Adaptive  stretching: 

2 min  | e , r/2 | 

2 max  |e,  r/2,  0| 


as 


■I 


for  < r/2 
for  > r/2 


(4) 


(5) 


Each  is  based  on  a 3 x 3 array  of  picture  elements  which  are 
illustrated  in  Figure  2. 
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3x3  Array 


a 

b 

c 

d 

e 

f 

g 

h 

i 

Figure  2.  Kernel  of  Pixels  Used  in  the  Calculations, 

Illustrating  the  Notations  Used  in  Equations  1 
through  5 . 

The  first  circuit.  Figure  la,  performs  the  Sobel 
operator  detecting  edges  in  two  dimensions.  The  processor 
architecture  is  arranged  in  the  form  of  a two-dimensional 
transversal  filter  with  impulse  response  for  the  two  edge 
components . 
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(6b) 


Using  these  two  components,  both  the  absolute  magnitude  of 
the  operator,  eq. (1) , and  the  edge  direction,  tan  0 = S /S 

x y 

are  directly  available.  The  effectiveness  of  the  weighting 
techniques  is  shown  in  Figure  3.  The  two  edge  components  Sx 
and  Sy  then  are  applied  directly  to  a CCD  absolute  magnitude 
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Sy  - 1/4,  0,  1/4 
(not  to  same  scale) 


V -v 


* * r 


W—  X 


Sy  = -1/8,  -1/4,  -1/8 


Figure  3.  Impulse  Response  of  the  2-D  Filter 


operator  and  a charge  summer. 


i 


The  performance  of  the  CCD  edge  detector  is  illustrated 
in  Figure  4 where  an  original  black  and  white  test  pattern 
(4a) , the  computer  simulated  Sobel  (4b) , and  the  output  of 
the  CCD  processor  (3c)  can  be  compared.  The  clock  rate  for 
this  demonstration  was  4 kHz,  limited  primarily  by  the  test 
facilities.  For  comparison,  the  CCD  Sobel  operation  of 
other  optical  images  is  given  in  Figures  5-7.  Our 
evaluations  [2]  indicate  that  at  these  clock  rates,  the 
operation  has  an  accuracy  and  dynamic  range  equivalent  to 
four  bits.  We  are  currently  unable  to  examine  a larger  gray 
scale  due  to  the  access  time  of  the  processed  data  from  the 
commercial  refresh  memory  we  are  using. 


We  have  spent  considerable  effort  developing  a real 
time  processing  capability  to  operate  the  CCD  processor  from 
a commercial  vidicon  camera,  the  Cohu  Model  No.  7120  13]. 
The  basic  data  rate  required  for  this  is  approximately  7.5 
MHz.  We  are  currently  operating  our  CCD  processor  at  2 MHz, 
which  results  in  a slightly  unsymmetr ical  Sobel  operation  as 
shown  in  Figure  8.  The  frame  rate  is  equivalent  to  60 
fields/sec  with  512  lines  as  in  standard  television,  however 
the  pixel  resolution  in  the  horizontal  direction  is  degraded 
by  approximately  a factor  of  3.  We  have  tested  the  circuits 
at  these  rates  with  a variety  of  images  and  our  intention  is 
to  increase  the  effective  data  rate  in  the  next  phase  of  the 
program  to  achieve  truly  symmetrical  operation. 


The  edge  detection  circuit  described  above  is  basically 
an  important  demonstration  of  two-dimensional  nonlinear 
processing.  Our  second  test  chip  which  performs  the 
operations  described  in  eqs.(l)  through  (5)  is  aimed  at 
demonstrating  adaptive  functions  based  on  the  local  mean  or 
average.  As  such,  the  prime  operators  are  the  edge 
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(b)  Computer  simulation 


(c)  Output  for  CCD  processor 


Figure  4.  Example  of  Sobel  Edge  Detection  Using  CCD  Processor. 
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(a)  Original  image 


(b)  Ce'jputer  simulation  (c)  Output  for  CCD  processor 


Figure  5.  Example  of  Sobel  Edge  Detection  Using  CCD  Processor. 
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Figure  6.  Example  of  Sobel  Edge  Detection  Using  CCD  Processor. 


Figure  8.  Example  of  the  Operation  of  the  CCD 

Sobel  Processor  Operating  in  Real-Time 
from  a Commercial  Vidicon. 
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detection,  local  averaging  and  the  delayed  original  images. 
Each  of  the  other  algorithms,  the  unsharp  masking,  the 
binarizer  and  the  adaptive  stretch  are  arithmetic 
combinations  of  these.  The  original  image,  the  Sobel  and 
the  3x3  mean  derived  from  the  second  chip  are  illustrated 
in  Figure  9 for  a regular  test  pattern.  Examples  of  the 
operation  on  a true  optical  image  are  shown  in  Figure  9. 
Each  function  described  in  eqs.(l)  through  (6)  (and  included 
in  Test  Chip  II)  has  been  tested  and  we  estimate  the  overall 
performance  to  be  equivalent  to  approximately  4 bits. 
Testing  of  linear  combinations  of  the  operators  described  in 
eqs.(4)  through  (6)  has  not  been  completed  at  the  full  video 
rates,  and  this  effort  is  currently  proceeding.  We 
anticipate  no  significant  problems  in  this  area. 

New  Concept  Development 

In  addition  to  the  above  work,  we  have  started  concept 
development  and  analysis  of  a third  test  chip  to  perform 
statistics,  including  a 7x7  median  filter,  an  analog 
histogrammer  (including  a mode  and  standard  deviation 
filter,  a 7 x 7 programmable  processor,  and  a number  of 
bipolar  fixed  filters) . This  work  will  continue  into  the 
next  phase  of  the  program  when  the  detailed  design, 
simulation  and  initial  processing  will  be  undertaken. 

Development  of  a Real  Time  Demonstration  Unit 

As  part  of  our  effort  to  interface  the  currently 
developed  processors  with  a commercial  video  camera,  we  are 
pursuing  the  development  of  a small  real  time  demonstration 
unit  which  will  include  the  necessary  analog  CCD  delays,  the 
clocks  and  drivers  for  our  processor,  the  CCD  processors 
themselves  and  a small  video  display  unit.  This  work  is 
well  underway,  most  of  the  interface  circuitry  having  been 
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Figure  9.  Illustration  of  Performance  of  Test  Chip  II. 


designed,  and  we  plan  to  have  the  complete  unit  available  in 
the  next  phase. 

Conclusions 


During  the  previous  phase  of  this  program,  we  have 
developed  CCD  integrated  circuit  processors  which  perform 
two-dimensional,  nonlinear  and  adaptive  operations  at  speeds 
in  excess  of  two  orders  of  magnitude  higher  than  general 
purpose  computers.  Our  evaluations  of  this  circuit  to  date 
indicate  that  it  will  perform  as  predicted  [4]  and  can  be 
interfaced  directly  to  the  optical  sensors,  and  as  such  lead 
directly  to  the  development  of  truly  smart  sensors. 
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4.2  Statement  of  Work  For  Follow  On  CCD  Circuitry 
Harry  C.  Andrews 


Introduction 


This  statement  of  work  is  designed  to  outline  various 
research  and  development  tasks  to  be  performed  by  the  Hughes 
Research  Laboratories  in  cooperation  with  the  USC  Image 
Processing  Institute  under  funding  from  the  Defense  Advanced 
Research  Projects  Agency  (DARPA) . It  is  anticipated  that 
there  will  be  four  projects  spread  over  the  two  years  of  the 
contract  (October  1,  1977-September  30,  1979)  ranging  from 
hardware  demonstrations  to  research  studies.  The  four  tasks 
are  delineated  below. 

Real  Time  TV  Demonstration 

At  the  request  of  DARPA,  task  1 will  be  to  develop  a 
portable  self  contained  television  demonstration  of  the  CCD 
circuits  developed  on  this  and  past  efforts  for  USCIPI  and 
DARPA.  Such  circuits  could  include  the  Sobel , Circuit  II, 
and  possibly  the  7x7  circuits  to  be  concurrently  developed 
over  this  current  two  year  effort.  It  is  anticipated  that 
the  first  real  time  TV  demonstration  unit  be  available, 
nominally  6 months  after  initiation  of  this  contract. 
However  in  the  interest  of  simplicity,  real  time  TV  can  be 
interpreted  to  mean  processing  every  other  field  at  regular 
TV  rates  thereby  implementing  the  chip  operation  on  nominal 
262  x 525  image  planes.  However  it  is  recommended  that  no 
scan  rates  be  changed  in  order  that  flicker  and  other 
artifacts  not  be  developed. 
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Completion  Testing  of  Circuit  II 


Circuit  II,  developed  on  previous  efforts,  should  be 
tested  as  the  second  task  of  this  contract.  This  testing 
should  result  in  5 kHz  demonstrations  of  the  various 
functions  designed  into  the  circuit.  In  addition  the 
question  of  SNR  and  dynamic  range  should  be  answered. 
Attempts  to  achieve  1 part  in  256  dynamic  range  accuracy 
should  be  made.  It  is  suggested  that  even  if  1 part  in  256 
can  be  achieved  on  a nonlinear  but  monotonic  scale,  this 
will  be  a significant  advance.  Nonlinear  computer  generated 
test  patterns  are  suggested  for  this  experiment. 

7x7  Circuits 

This  task  will  probably  require  the  greatest  amount  of 
CCD  innovative  design  capability  and  probably  represents  the 
heart  of  the  follow-on  effort.  A series  of  circuits,  each 
7x7  pixels  in  extent,  is  to  be  designed  fabricated  and 
tested  (certainly  at  5 kHz  and  hopefully  at  TV  rates) . A 
variety  of  7 x 7 functions  are  to  be  implemented  to  achieve 
a)  fixed  bipolar  linear  filtering,  b)  operator  programmable 
bipolar  linear  filtering,  c)  data  programmable  linear 
filtering,  and  d)  texture  related  7x7  parameter 
development . 

a)  The  fixed  bipolar  linear  filtering  circuit  will  be 
designed  to  test  both  the  ability  to  implement  bipolar 
coefficient  processes  and  to  measure  the  accuracy  of  fixed 
coefficient  numeric  values.  The  filter  suggested  will  be 
biquadrant  symmetric  defined  by  16  degrees  of  freedoms  as: 
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b) 

The  operator 

programmable  bipolar 

linear  filter  is 

designed  to  test  the  capability  of  switching  filters  in  the 
middle  of  a processing  mission.  Electronic  signals  are 
anticipated  for  use  in  this  switching  process.  The  above 
filter  could  be  one  of  the  circuits  while  another  circuit  of 
interest  might  be  another  biquadrant  symmetric  filter 
defined  by: 
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Let  a be  a controllable  parameter  (set  to  unity  for  this 
circuit) . 

c)  The  data  programmable  linear  filter  is  one  which 
changes  its  structure  as  a function  of  the  current  data 
within  the  7x7  window  under  processing.  The  observant 
reader  will  note  that  the  above  filter  is  unipolar  (except 
for  the  center  pixel) . If  the  unipolar  portion  is 
implemented  separately,  then  the  size  of  a can  be  controlled 
by  the  unipolar  (low  pass)  portion  of  the  filter.  Letting  a 
be  a monotonic  function  of  this  low  pass  signal  then  will 
cause  greater  unsharp  masking  (edge  enhancement)  in  the 
brighter  regions  and  less  edge  enhancement  in  the  less 
bright  regions. 
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d)  The  texture  related  7x7  circuit  need  not  have  the 
precision  in  coefficient  definition  given  by  the  above 
filters.  However,  what  the  texture  circuit  lacks  in 
coefficient  accuracy,  it  more  than  makes  up  for  in 
cleverness,  variety  and  nonlinear  on  chip  decisions.  In 
order  that  as  much  creativity  and  unique  CCD  design  be  made 


available  for  this 

texture  circuit,  a variety  of 

suggested 

operations  are  listed  below  for  design  study 

fabr ication . 

prior  to 

1) 

median  filter 

2) 

mode  filter 

3) 

dispersion  filter 

4) 

standard  deviation  filter 

5) 

checkerboard  textures 

The  median  filter  outputs  the  median  grey  level  of  the 
49  pixels  within  the  7x7  window. 

The  mode  filter  outputs  the  mode  grey  level  of  the  49 
pixels  within  the  7x7  window. 

The  dispersion  filter  outputs  the  height  (amplitude  or 
frequency  of  occurance)  of  the  mode  of  the  49  pixels  within 
the  7x7  window. 

The  standard  deviation  filter  outputs  the  square  root 
of  the  variance  within  the  7x7  window.  (consider  sum  of 
absolute  differences  rather  than  squares  or  square  roots) . 

Checkerboard  textures  output  the  value  of  the 
correlation  of  the  checkerboard  with  the  data  within  the 
7x7  window.  There  are  8 such  checkerboards  to  be 
considered.  See  Figure  1. 
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It  is  hoped  that  after  careful  study  of  these  5 
categories  of  texture  circuits,  a clever  implementation 
technique  will  emerge.  Certainly  a large  subset  of  these 
suggestions  should  be  implemented  in  the  7x7  texture  chip. 

Seqmentor  Study 

The  final  task  for  this  two  year  effort  should  be  the 
study  of  possible  implementation  of  a focal  plane  (CCD) 
image  segmentor . This  should  include  the  combination  of 
outputs  of  previous  circuits  and  the  definition  of  areas  in 
need  of  further  technological  development.  This  task  will 
be  a continually  evolving  process  over  the  two  years  and 
should  include  dialogue  between  USCIPI  and  HRL  personnel. 
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5.  Hardware  Activities 


This  section  of  the  semi-annual  report  is  devoted  to 

/ 

recent  Institute  hardware  activities  dealing  with  digital 
imagery  equipment. 


5.1  Hardcopy  Acquisition 
Harry  C.  Andrews 


The  Image  Processing  Institute  has  recently  acquired  a 
Dicomed  D47  color  and  monochrome  hardcopy  device.  This 
equipment  now  allows  us  to  output  high  resolution  digital 
imagery  (4096  x 4096) : 


a) 

4" 

x 5" 

Polaroid 

b) 

4" 

x 5" 

sheet  film 

c) 

35 

mm . 

film  strips 

d) 

35mm. 

cine  film 

The  device  has  already  proven  quite  useful  and  has  been 
calibrated  with  our  own  digital  test  chart.  Figure  1 
presents  this  test  chart  at  varying  resolution  display 
parameters.  Use  of  the  display  includes  traditional 
photographic  work  as  well  as  digitally  generated  halftone 
screens  for  subsequent  processing  in  our  hybrid 
optical-digital  laboratory. 
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(a)  lx  and  2x  magnification 


(b)  4x  magnification 


Figure  1.  Hardcopy  Display  of  USCIPI  Test  Chart. 
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5.2  The  RTTV  at  ARPA 
Harry  C.  Andrews 


The  Image  Processing  Institute  has  installed  a real 
time  television  (RTTV)  display  at  ARPA  headquarters  in 
Washington,  D.C.  The  device  is  capable  of  storing  and 
refreshing  a 256  x 256  x 8 bit  monochrome  image  and  has  an 
8x8  brightness  function  memory.  The  software  which 
supports  the  system  is  fairly  simple  and  allows  ARPA 
personnel  the  ability  to  access  remote  image  files  for 
up-to-date  results  of  various  research  projects.  Figure  1 
presents  a picture  of  the  device  driving  a Sony  monitor  and 
our  Advent  display. 

It  is  expected  that  CMU  will  soon  have  the  software  for 
access  by  ARPA  personnel  to  their  image  files.  Other  ARPA 
contractors  wishing  to  avail  themselves  of  this  novel  use  of 
the  ARPANET  are  invited  to  pursue  the  requirements  with 
their  contract  monitors  and  with  USCIPI  personnel. 
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Figure  1.  The  RTTV  Display  (before  shipment  to  ARPA 
in  Washington,  D.C.) 
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6 . Recent  Ph  ,D.  Dissertations 

This  section  includes  all  dissertations  completed  since 
the  last  reporting  period.  The  one  listed  here  reflects  the 
results  of  applying  the  mathematical  model  of  the 
psychophysics  of  vision  to  image  compression.  Results  are 
extremely  encouraging  and  experimental  demonstrations  of 
24:1  compressions  without  noticeable  degradations  are 
obtainable.  The  rate  distortion  analysis  also  indicates 
that  this  model  produces  results  close  to  the  theoretical 
limit. 

Detailed  results  of  this  research  topic  are  described 
in  detail  in  USCIPI  Report  #790. 


6.1  Digital  Color  Image  Compression  in  a Perceptual  Space 
(Supported  by  WPAFB  under  Contract  F-33615-77-C-1016) 

Charles  F.  Hall 


Simple  mathematical  models  are  developed  from  the 
physiological  and  psychophysical  traits  of  the  human  visual 
system.  Expressions  for  the  statistical  characterization  of 
these  models  are  obtained.  When  used  as  a preprocessor,  the 
models  are  shown  to  produce  images  which  are  statistically 
compatible  with  the  underlying  assumptions  necessary  to 
solve  the  parametric  rate-distortion  equations.  The  derived 
power  spectrum  equations  were  used  to  code  black  and  white 
and  color  images  with  a quality  superior  to  previous 
results.  In  addition,  it  is  shown  that  the  preprocessor 
, r 1 . es  a "perceptual  space"  in  which  normalized  mean 
. • r ts  an  effective  image  quality  measure. 
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